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Chapter 1 


BACKGROUND: TURBULENT FLOW COMPUTATION METHODS 


1 . Methods of Computing Turbulent Flowa ; Clasalfication 

A few years ago, the author and two of his colleagues wrote a paper which 
attempted to classify methods of dealing with turbulent flows (Kline et al. 
(1978)). This paper is reviewed and extended here as a means of setting the 
main subject of this report in context. 

There are two sub-areas that need to be dealt with in classifying methods 
of computing turbulent flows. These are the method by which the fluctuations 
are treated and the manner in which the geometry of the flow is handled. 
These arc, of course, coupled to some extent, but it is useful to separate 
them for purposes of this work. We shall take up the problem of dealing with 
the turbulence first. According to the classification scheme in the paper 
cited above, there are five broad classes of methods of dealing with the 
turbulence; there are also subclasses of each. The fivf> major categories are: 

i) Correlations . These are the familiar Correlations that give the 
nondlmeuBional skin-friction coefficient as a function of the Reynolds number, 
Nusselt number as a function of Reynolds and Prandtl numbers, etc. They are 
extremely useful, but very limited. Their applicability is especially limited 
in high-technology applications in which the geometry plays an important role 
in the fluid dynamics (such as airfoils); for such problems, a new set of 
correlations would be needed each time the geometry of the device is changed. 

ii) Integral Methods . In these methods the equations governing the 
fluid dynamics (which may be the equations used on level (Hi) below) are 
integrated over at least one coordinate direction. This decreases the number 
of independent variables and greatly simplifies the mathematical problem to be 
solved. These methods allow considerable use of experimental deta and 
physical insight and have proven quite useful. One of their principal 
drawbacks is that they need to be reworked when a new type of flow is to be 
computed . 
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til) Reynolds-Averaged Equations . In this approach, one averages the 
Navier-Stokes equations over either time, homogeneous directions in the flow, 
ov an ensemble of essentially equivalent flows . When averaging of any of 
these kinds is performed, the equations describing the mean field contain 
averages of products of fluctuating velocities, and there are fewer equations 
than unknowns — the well-known closure problem. In fact, the set of equations 
can never be closed by further averaging; a closure assumption or, what is the 
same thing, a turbulence model has to be introduced. The closure assumption 
must represent the unknown higher-order average quantities in terms of the 
lower-order quantities that are computed ejcplicltly. This subject is 
undergoing a rapid expansicn at the present time. It is also likely thrt this 
level should be broken into sublevels or separate levels. 

iv) La rge Eddy Si mulation . In this approach, the equations are averaged 
over a small spatial region. The object is to remove the ^rnall eddies from 
the flow field so that an equation for the large eddies is derived. The ef- 
fects of the small eddies on the large ones is then modeled . This is one of 
the principal subjects of this report and is discussed in considerable detail 
below. 

v) Full Simulation . This is the numerical solution of the exact 
Navier-Stokes equations. The only errors made are numerical ones which, with 
care, can be kept as small as desired. By its nature, this approach is 
limited to low Reynolds numbers. This is the other principal subject of this 
report and will be covered in detail below. 

Currently, computations at levels (Iv) and (v) are limited to people with 
access co very large, fast computers. They are not suitable for engineering 
design at present and we anticipate that it will be some time before they will 
be (if ever). We call levels (iv) and (v) together higher-level methods of 
turbulence computation— hence the title of thlr* report, 

A significant point about this classification scheme is that calculations 
on any levels can be used to generate information that can be used on the 
lower levels. In applications, engineers commonly use methods at level (ii) 
or (ill) to I roduce correlations from which the design is actually done. 
Large eddy simulation (LES) can be used to produce information chat can be 
used in modeling for Reynolds-averaged calculations . LES could be used in 
principle at the lower levels as well, but there is little need for this 
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application* Full simulation can be used Co test models for both the 
Reynolds-averaged equations and large eddy simulation. this will receive 
considerable attention in this report. 

It should be noted that the nomenclature we have used for classifying 
methods differs from chat of Schumann et al. (1980). What we have called 
higher- level simulations they called direct 'Simulation, and they did not make 
the distinction between levels (Iv) and (v). We believe the discincCion 
important and prefer the nomenclature used in this, report. 

The second type of classification of methods of computing turbulent flows 
concerns the treatment of Che geometry. This scheme contains just two cate- 
gories : 

a) Full Field Methods . In this approach, the same set of equations is 
applied everywhere in the flow field. This has the great advantage of not 
requiring any kind of matching In the Interior of the flow and of being easier 
to program for computer solution, The principal drawback is that fine meshes 
are needed in some regions of the flow (such as near Che boundaries and in 
shocks), and this can make the cost very high. 

b) Zonal Methods . In zonal methods the flow is considered as a 
collection of "modules ," and each module or zone Is treated by a separate 
method. The most common example of this kind of method is the division of 
flows over bodies into boundary layers and potential flows which are treated 
by separate methods. In zonal methods, Che solutions in the various zones 
have Co be matched at their common boundaries, by an iterative process Chat 
may or may not converge well. The modules can be Created by different 
methods. Thus one can use an integral method for Che boundary layer and the 
full partial differential equations for the outer flow. 

The classification scheme given here differs a little from the earlier 
one of Kline et al. (1978). We believe that the current scheme represents an 
Improvement in clarity. We have found it useful, and it will be one of the 
ways in which various methods will be compared in the 1980-81 SCanford-AFOSK 
Symposium on the Computation of Complex Turbulent Flows. 
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2 . 


Classification of Turbulent Flows 


An issue that is quite separate from that of how turbulent flows are 
computed is that of trying to classify the flows. In a field as complex as 
this, any classification scheiae is inexact, but it is better than having no 
scheise at all. Thus we shall classify flows according to the phenomena that 
occur in them. This scheme is not new and contains three categories: 

a) Homogeneous Flows . In th!!se flows the state of the fluid is the 
/jame at every point in space; they develop in time. There is a limited number 
of flows of this kind; the experimental data for them have been reviewed re- 
cently by the author (Ferziger (1980)). In homogeneous flows without mean 
strain or shear, the turbulence decays with time; when mean strain or shear 
are applied, the kinetic energy of the turbulence may increase with time. The 
mechanism by which the turbulenee length scales increase in these flows is not 
well vinderstood. 

b) Free Shear Flows . It is well known that free shear flows are 
extremely unstable. The laminar mixing layer is unstable with respect to 
disturbances over a wide range of wavelengths. The instability is of the 
Kelvin- Helmholtz type in which the perturbation grows rapidly. There is 
controversy about the precise mechanism of growth of the turbulent free shear 
layer, but it seei^s clear that there are large coherent regions of 
concentrated vorticlty in all of these flows. The concentrations of vorticity 
cause strong large-scale mo^4ons within the flow and the vorticlty tends to 
agglomerate further. The controversy centers on the nature of the 
agglomeration, c,f. Roshko (1978) and Chandrsuda et al. (1977). 

A subclassification of these flows is necessary. In the mixing layer 
(the simplest type of free shear layer), the velocity difference across the 
layer remains fixed as the layer develops. As a result, the layer grows lin- 
early in space or time, InJefinitely In other flows, for example, jets and 
wakes, Che velocity differences are reduced as Che flow develops and the 
turbulence tends Co weaken in the downstream direction. 
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e) Wall-Bounded Flows . The effect of a wall on a shear layer Is to 
prevent (or at least reduce) the large-scale motions described in the previous 
paragraph and thus inhibit the shear layer from growing so rapidly. Thus, 
boundary layers and related flows grow less rapidly aa*l have lower turbulence 
levels than do free shear layers . Another , weaker , mechanism of turbulence 
production takes over. This mechanism is less well understood than that of 
the free shear layer and, perhaps for that reason, seems much more complica- 
ted. It is known to involve Che presence of thin regions of high- and low- 
speed fluid that exist close to the wall and which are long in streamwlse 
extent (Kunstadler et al. (1967), Kim (1969)) and large-scale motions of the 
outer part of the boundary layer, but several details remain to be filled in. 

A further extension of this classification scheme was given by Bradshaw. 
His view is that the mean turbulent flows can be thought of as a combination 
of “normal" strains — the mean ettains that occur in the “standard” flows-' 'and 
“extra" rates of strain. There are many extra rates of strain. Some of them 
are: curvature, rotation, lateral divergence (in axisymmetric flows), 
buoyancy, blowing or suction, and wall roughness. Although these effects 
generally appear as small terms in the equations, they have profound effects 
on the structure of the turbulence and. Indirectly, on the behavior of the 
flow as a whole. Therefore, they are very important, and we shall devote part 
of this report to investigating their effects on tiarbulent flows. 

Finally, it should be noted that some complex flows may be of one type in 
one region and another type In another region. In particular, in flows with 
separation, wall boundary layers may become free shear layers and vice versa. 

3 . A Short History 

There are no known analytical solutions of the Navler-Stokes equations 
for turbulent flows, and it is unlikely that there ever will be any. This 
fact, plus the obvious technological importance of turbulent flows, is the 
reason for the development of computational methods of predicting turbulent 
flows . 

Prior to 1960, computers had too little capacity to do anything more than 
solve the ordinary differential equations of Integral methods or the partial 
differential equations for simple, two-dimensional potential flows. Progress 
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in this period was larp,eXy i'escricted t*o the computerization oS methods that 
had been carried out on der^k calculatora up to that time. 

At< computer's grew in Bophiaticatlon, eo did the problems for which people 
sought solutions. The 19608 saw the development of good boundary layer meth- 
ods based the use of both Integral methods and partial differential equations 
levels (li and ill of the above scheme). The 1966 Stanford Conference (Kline 
(1968)) marked a milestone in this development. At that time, people were 
beginning to solve the Reynolds-averaged Navier-Stokes equations using simple 
models for relatively simple flows. Through the 1970s, the sophistication of 
the models grew, as did the complexity of the flows that researchers were 
willing to try to compute. 

The first applications of what we have defined as higher-level methods 
Were made by the meteorologists. That field has needed models for predicting 
the world’s weather patterns for a long time* As soon as computers were large 
enough, meteorologists tcied global weather simulations. The first three- 
dimensional attempt at this of which the author is aware is Chat of Smago- 
rinsky (1963); this paper presented a model that will be used extensively 
later in this report. Ihe grid systems used in these early calculations were 
necessarily very coarse, and the method used was necessarily what we have 
called large eddy simulation. Improvements in computers have allowed the use 
of finer grids, but the grids are still coarse compared to what is desired, 
this situation, unfortunately, will not change in the foreseeable future. 
Hence, subgrid-scale modeling will remain an important issue in meteorology 
(and oceanography) for quite some time. 

Tlie first computation of a flow of engineering interest was the simula- 
tion of channel flow by Deardorff, a meteorologist, in 1970. In this landmark 
paper, he laid out many of the foundations of the field. Improvements in his 
methods were made by Schumann (1973) and Grotzbach (1976). Tlie latter and 
their group at Karlsruhe have subsequently extended the method to the 
computation of annular flows, the inclusion of heat transfer, and the 
inclusion of the effects of buoyancy. 

The author's group at Stanford, which is jointly led by W. C. Reynolds, 
began work in higher-level simulations in 1972. Their objective was to put a 
sound foundation under the method of large eddy simulation by computing simple 
flows first. It was felt that in this way the fundamentals of the subject 
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could slowly be put In ordec. Die first flows chosen for study were the 
honogeneous turbulent flows, and quite Vot was learned about nunerlcal 
methods and subgrid-scale modeling (Kwak et nl, (1975) and Shaanan et al. 
(1975)), When the group felt that the techniques for tlie simulation of homo- 
geneous flows were well developed, it was decided to go on to the study of 
flows which are inhomogeneous in one coordinate direction. Die simplest such 
flows are the mlKing layer and channel. Die fully developed mixing layer was 
computed by Mansour et al. (197B), transitic'< in the mixing layer was studied 
by Cain et al. (1981), and the channel flow was studied by Moin et al. (1978) 
and Kim and Moln (1980,1981). 

Almost from the beginning it was realized that the effort in computing 
flows would have to be accompanied by an effort at developing better models 
for creating the small scales (subgrid scale models) or at least understanding 
the models chat are in ‘.isc. The method of using direct simulations for this 
purpose was developed by Clark et al. (197b) and extended by McMillan and 
Perzlger (1978), McMil’.an et al. (1980), and Bardina et al. (1980). 

It is clear that large eddy simulation will not be a method of direct 
engineering applicability for some time. For that reason, the major impact 
the method will have is in the improvement of the understanding of the physics 
of turbulent flows and in helping to develop, test, evaluate, and improve 
model's that are used in Keynolds-averaged methods. Recently, exact Simula-^ 
lions of compressible homogeneous turbulent shear flows and homogeneous 
turbulent shear flow with a passive scalar weie made in order to evaluate 
these models; cf Feiereisen et.al. (iS'81, and Shlrani et.al. (1981). 

A group under Leslie in London has been active in the field since 1975. 
Their early work centered on the understanding of subgrid scale models (Love 
and Leslie (197b) and Leslie and Quarini (1979)). Since then they have cimu- 
iated the mixing of a passive scalar in homogeneous isotropic turbulence 
(Antonopoulos (1981)). 

A number of French groups have studied subgrid scale models from a theo- 
retical point of view and have made several contributions in this area. 

Orszag and coworkers have been working since 1970 on the direct simula- 
tion of turbulent flows. Their early work centered on the prediction of 
homogeneous Isotropic turbulence (Orszag and Patterson (197i)), and more 
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recencly they have become interested in the study of transition in wall- 
bounded flows (Kells and Orszag (1979), Orszag and Patera (I960)). The main 
interest of this rroup has been in the development of numerical methods (they 
are responsible for the widespread use of spectral methods in this field), on 
the study of turbulence theories, and on the prediction of transition* 

Kiicy and Metcalf (19BU) have made direct simulations of free shear 
flows, ‘ineir efforts have been directed at the simulation of fully developed 
Wakes at relatively low Reynolds numbers, which may be thought of as the last 
stages of the decay of a turbulent wake. 

Rogaiio (197B, 19B1) has made extensive direct simulations of all of the 
homogeneous turbulent flows. Ills results are an important resource for 

modelers . 

4 . Outline of Tliis Report 

in Chapter II, we shall consider the fundamentals of large eddy 
simulation and compare the various approaches to it. 

In Chapter III we shall discuss the subgrid scale models required by 
large eddy simulation. We shall also study the use of large eddy simulation 
in the development of models for the Reynolds-averaged equations and the 
application of full siraulation to the testing of both subgrid scale and 
Reynolds-averaged models. 

In Chapter IV we shall discuss the numei'ical methods used in large eddy 
and full simulation. Since the numerical methods used are almost always 
somewhat tailored to a particular flow, we shall just touch on some of the 
special-purpose methods in this chapter. the latter methods will be 
considered in more, detail in the chapters In which the flows are described. 

Chapter V will be devoted to the discussion of the simulation of 
homogeneous flows. Tlie flows will be categorized, and the numerical methods 
needed for some of the cases will be described, along with physical 
descriptions of the flows. We shall give the results from both full and large 
eddy simulations of these flows and show how they can be applied to the 
testing and development of models. This chapter contains a considerable 
amount of recent work. 
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Free shear flows will be considered in Chapter VI. Hie bulk of the chap- 
ter will be devoted to the mixing layer, which has been the principal focus of 
attention in this area, but we shall also look at wake simulations. 

Chapter VIX will be concerned with wall-bounded flows. Most ot the at- 
tention will be given to channel flow, but some discussion of recent work on 
the boundary layer will also be given. Particular attention will be given to 
terms which have not been measured in the laboratory. 

In Chapter Vlll we shall briefly cover applications of large eddy simula- 
tion and full simulation that have not been given in the previous chapters. 
Hie most important of these applications are in muteoroiogical and other 
environmental flows. However, a few applications have been made to other 
laboratory flows, and these will be briefly covered as well. 

T»ie concluding chapter, IX, will discuss some directions in which the 
work is proceeding and what can be expected from higher-level simulations in 
Che next few years. 

Tills report will give greater emphasis to work done in the author's group 
than to that of other groups. Tlie reader is reminded that this is a conse 
quence of greater familiarity with his own work and that of his colleagues and 
is in no way intended to imply that work done elsewhere is any less important. 
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Chapter II 

FOUNDATIONS OF LARGE EDDY SIMULATION 


1 B Rationale 

It is generally believed that the largest eddies dominate the physics of 
any turbulent flow. The differences between the large and small eddies can be 
summarized as follows: 

a) The large eddies Interact strongly with the mean flow. The small 
eddies are created mainly by nonlinear interactions among the large eddies. 

b) Most of the transport of mass, momentum, energy, and (in flows 
containing more than one species) concentration is due to the large eddies. 
The small eddies dissipate fluctuations of these quantities but affect the 
mean properties only slightly. 

c) The structure of large eddies is very strongly dependent on the geom- 
etry and nature of the flow. They are usually vortical, but their shapes and 
strengths are flow dependent. The small eddies are, on the. other hand, much 
more universal. 

d) Due to their dependence on the geometry, the large eddies are highly 
anisotropic, llie small eddies are much more nearly isotropic and, therefore, 
universal. 

e) The time scales of the large eddies approximate the time scales of 
the mean flow. For flow past a body, the large eddy scale is approximately 
the dimension of the body divided by the free stream velocity. The small 
eddies seem to be created and destroyed much more quickly. 

An important consequence of these properties is that the large eddies 
should be harder to model than the small ones. Also, as they vary so much 
from flow to flow, one should not expect to find a model for the large eddies 
to be universal. There is hope, however, that one might be able to find a 
useful model of the small eddies. 

This leads to the concept of large eddy simulation. In this approach, 
the large structures in the flow are computed explicitly and the small ones 
are modeled. This method should have advantages over methods in which all of 
the eddies are modeled. 
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TliCHf argumonts provide tUe raLlonuie lor large eddy aimulatlou. 
However, not all of tlic premises given above hold in ail Hows. Tliey seem to 
hold iJi Iwmogeneous turbulent flows and in free shear flows, in wall-bounded 
Hows, liowever, the structures responsible for much of the luomentum transport 
(and, presumably, tiie transport of the other properties as well) may be quite 
small, especially in the region close to the solid boundary. Speci»f cwire is 
necessary in these flows, this will be discussed further in Chapter Vll. 

If, for now, we accept tlie statements made above as correct, it follows 
(hat large eddy simulation ought to have a number of advantages over Reynolds 
or time average metliods. The most significant advantage is that much of the 
actual transport of mass, momentum, energy, and species is compvited 
explicitly, and the portions of these fluxes which need to be, modeled are tiuich 
smaller than what is modeled in the Reynolds-averaged equations. 
Consequently, tl>e overall results are less sensitive to modeling inaccuracy in 
large eddy simulation than in other approaches. Tlie probability of iinding a 
widely appilcalde model should be much higher. 

n»e principal disadvantage of large eddy simulation relative to Reynolds- 
averaged methods Is that the computations are necessarily three-dimensional 
and time-dependent. 'nils means that the cost is much higher. In fact, the 
cost is currently high enough that, except for the simplest Hows, use oi the 
method Is restricted to people with access to large amounts of time on very 
la rge computers . 

ilte first tiling that one needs to do In developing large eddy simulation 
is to define the large scale component of the flow fields the portion which is 
to be computed expiicitly. Ihere are two common approaches to doing this, 
tliey will be described and compared in the next two sections. Tlie, remainder 
of the chapter will present the UiS equations and describe the puramett 
tradeoffs that must be faced in large eddy simulation. 

Tlie equations for the large scale field alw,jys contain terms which 
involve the small scale field, which Is not computed. These terms play the 
same role in the large eddy equations as the Reynolds stresses play in the 
Reynolds-averaged equations. They are therefore called tlie subgrid scale 
(SUb) Reynolds stresses, and they must be modeled. A discussion of subgrid 
scale models and a comparison of them with Reynolds-averaged turbulence models 
Is given In the next chapter. 
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2. Fllterlnij 


The Tirsc task in large eddy simulation is that of defining the large 
scale component of the flow field — tlie portion which the method will attempt 
;.o calculate. Tliere arc several ways of doing this mathematically. All are 
essentially equivalent to averaging the equations over a vHmall region of space 
or low-pass filtering the equations in Fourier space. llie starting point Is 
the incompressible Navier-Stokes equations: 


3u. » I ~ 
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( 2 . 1 ) 


which must be solved together with the continuity equation: 


du 

Jk 



( 2 . 2 ) 


For homogeneous flows, we prefer to define the large scale field (also 
culled the resolved field) by means of a convolution filter: 


in Fourier space, this has the form: 

u(k) “ G(k) u(_k) (2.4) 

Note tluit for this kind of filter C is a function only of the magnitude of 

Jc. 

A number of simple filters have been used. 111686 are illustrated in Fig. 

2.1. If the equations are simply integrated over a small control volume in 

space, we have the box filter i most finite-difference and finite-volume 
methods implicitly use this filtar. Its Fourier transform is also shown in 
the figure. Another common choice is a sharp cutoff filter in Fourier space, 
this is essentially the Fourler-space version of the first filter. Both of 
these filters have the difficulty that their Fourier transforms have negative 
regions, they also are difficult to differentiate. For this reason, we prefer 
to use the Gaussian filter. Its Fourier transform is also Gaussian, so it is 
well behaved in both configuration and Fourier space, and it can be differen- 
tiated as many times as one likes in both spaces. Tlie Gaussian is 


(2.5) 
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uc, in Kouricr space. 
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The iiumerlcal factors have been chosen to make the second moment of this flX- 

A 

ter the same as that of a box filter of width A. C(r) and U(k) are 
Fourier Inverses of each other when the variables are continuous, but not In 
the discrete case; In the latter case, a choice has to be made. tlie 
normullzatlon factor, A, has been left unspecified In Fq. (2.5) In order to 
admit the conservation property that the Integral of G(r) over all space be 
unity whether continuous or discrete quadrature Is used. 

Use of this type of filter was first suggested by Leonard (197'J), he 
showed how the concept could be generalized. It is sometimes useful to use 
expansions other than standard Fourier series. For example, Uhebychev 
polynomial expansions have been used (Orszag (197b), Kim and Moln (19b0)) as 
the basis for numerical methods. The filter can be defined In the space of 
the Index of the expansion functions; a sharp cutoff (Ignoring all components 
of the expansion beyond some specified N) Is the simplest possibility, but 
it is easy to construct Gaussian-like filters as well. 


When the filter (2.3) is applied to the Navler-Stokes equations (2.1) and 
the continuity equations (2.2), we have: 
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and 
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The difficulty comes from the nonlinear term. The approach taken by everyone 
In the field Is to write: 


ui - -F 

which causes the nonlinear terms to take the form: 


(2.9) 
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The first term Is entirely dependent on the large scale component of the field 
and is computable in LES. Hie small scale component of the velocity field, 
u|, is not computed, so the terms containing it need to be modeled, u^ is 
called the subgrid scale component of the velocity field, but this is a 
misnomer (in this formalism), because tive width of the filter (A) need not 
be related to the size of the grid on which the computations will be done, 
however, it has become standard nomenclature, and the set of terms involving 
the small scale velocity component. 


- uJ^Uj 4- u^u'j + ujij^ (2.11) 

are commonly called the subgrid scale (SCS) Reynolds stresses . Tliey must be 
modeled in large eddy simulations — hence the name s ubgrid modeling . We shall 
look at models for these terms in the next chapter. 

Tlie approach presented above is the cne favored by the author and his 
colleagues at Stanford. It decouples the definition of the large scale field 
from the numerical solution of equations that result. We favor this method, 
even though It is more cumbersome than the one given In the next section, be- 
cause we feel it provides more flexibility. This flexibility will be useful 
when we discuss methods of testing subgrid scale models in the next chapter. 

3. The Deardorf f-Schumann Approach 

An alternative to the method presented in the previous section is based 
on the recognition that we shall be solving the equations numerically. The 
computerprogram will be based on a set of discretized equations. It therefore 
makes sense to use an approach that arrives at the discretized equations as 
quickly as possible. The method originally presented by Deardorf f (1970) and 
extended by Schumann (1973) is one which accomplishes this. 

Hie Idea is to Introduce the grid on which the numerical computations 
will be done at the outset. Deardorff and Schumann used a staggered grid, 
which Is probably the best choice for solving the incompresCilble equations, 
but other grid systems could be used as well. The two-dimensional version of 
the staggered grid is shown in Fig. 2.2. One integrates each of the equations 


over an appropriate control volume; the control volume for the x-momentum 
equation la shown In Fig. 2.2. Hie resulting equations have the form (2.7) 
and (2.8), provided the operation represented by the overbar Is Interpreted as 
the volume average. Because the averaging operation Is defined relative to 
the grid, u^ Is defined only at the grid points. However, It Is convenient 
to think of Uj^ as constant within the control volume. This definition of 
the large scale velocity differs from the one presented In the previous sec- 
tion. Tlie two definitions are Illustrated In Fig. 2.3. 

The Oeardorff definitions lead to some convenient simplifications. In 
particular, one can assume that: 


UjU , “ u u (2.12) 

1 j 1 j 


and 


u{ - 0 (2.13) 

which are properties this approach shares with Reynolds-averaged modeling. 
Hie su'bgrld scale Reynolds stresses then, reduce to: 





(2.1A) 


Models are Introduced for and the discretized equations are similar co 
those commonly used on staggered grids. 

In Schumann's modification of this approach, the integrals of spatial 
derivatives are carried out analytically. This results in equations which 
contain Integrals over the surfaces of the control volumes. Hie difficulty 
with this approach is that four different types of averages appear: averages 
over the three types of faces of the grid volume and volume averages. These 
must be related In some way. Schumann Introduced several approximations that 
relate the surface averages to a single volume average, but the assumptions 
required are difficult to evaluate and may be questionable, especially at low 
Reynolds numbers. It Is not clear that this method has any significant advan- 
tages relative to Ueardorff's. 
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In the Deardorff-Schumann approach, the aubgrid scale velocity field la 
discontinuoua at the edges of the control volumes, and the behavior of the 
Bubgrid scale Reynolds stress as a function of position is not very smooth. 
Tills problem and the Increased flexibility in defining the filter are the 
primary reasons why we prefer the filtering approach to the one presented in 
this section. 


4. The Large Eddy Simulation Equations 

The equations of large eddy simulation are essentially (2.7) and (2.B). 
However, one needs to take into account Eqs. (2.10) and (2.11) as well. Also, 
one further modification is usually made. Tlie subgrid scale Reynolds stress, 
defined by Eq. (2.11), can be decomposed into the sum of a trace-free tensor 
and a diagonal tensor: 



(R 


ij 


T* 

1 


- T _ + fr 6^ .R. . 

ij J ij kk 


(2.15) 


Although the diagonal component of this tensor can be. modeled, there is no 
need to do so. When the decomposition (2.12) is substituted into the filtered 
Navler-Stokes equations (2.7), the diagonal component produces a term which is 
equivalent to the gradient of a scalar. It is similar to the pressure 
gradient term and can be combined with It. It is therefore advantageous to 
define a modified pressure: 


P 


P 


i\k 


(2.16) 


The filtered Navier-Stokes equations can then be written: 





3P 

7 ^ 


3^Ui 3 t^, 

j j j 


(2.17) 


Once a model for T^j has been Introduced, these equations are to be solved 
numerically along with the filtered continuity equation, which is repeated 
here for completeness: 


^ - 0 (2.18) 
j 
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5. Tradeoffs 


Xn any kind of flow computation, there are tradeoffs. Higher accuracy 
can always be by reducing the grid size and Increasing the number of mesh 
points. Tlie price Is paid In the form of Increased computer time. 

A similar tradeoff exists in large eddy simulation. Ideally, vt'e would 
like all the eddies In the large scale field to behave In the manner ascribed 
to large eddies at the beginning of this chapter and the small eddies to 
behave as they are supposed to. This separation of large and small eddies is 
possible only at high Reynolds. At sufficiently high Reynolds numbers^ the 
turbulent energy spectrum contains an Inertial subrange In which there .is 
essentially no turbulence production or viscous dissipation. The eddies which 
are larger than those In the subrange (l.e., lie at lower wavenumbers) behave 
like "large eddies", and those that lie at wavenumbers below the subrange are 
"small eddies." Since the width of the filter (A) Is supposed to mark the 
boundary between the two classes of eddies, the ideal is to choose the filter 
width such that the corresponding wavenumber (v/A) lies in the subrange. If 
this is the case, large eddy simulation should be successful. 

There are, of course, cUfflcultles that we need to address. Hie 
principal of these are: 

a) Tlie size of the physical domain considered in the calculation needs 
to be sufficiently large to hold the largest eddies. We also wish the filter 
size to be such that all of the "small" eddies lie in or below the subrange. 
Finally, the computational grid size must be smaller then the filter width 
(this Is discussed in Chapter 4). These requirements set the number of mesh 
points required In each coordinate direction. It is not unusual to find that 
the number of mesh points needed to meet these requirements Is much greater 
than the available computer resource will allow. We are then forced to use a 
filter width which lies outside the subrange. 

b) At low Reynolds numbers there is no subrange in the turbulence spec" 

trum. 

In either case, we are forced to use a filter width which does not lie 
within the Inertial subrange of the turbulence spectrum. It has been argued 
by some that one should not do this. We believe that it is reasonable to do 
large eddy simulation under these circumstances. However, the model may need 
to be changed to account for the fact that the cutoff is not in an inertial 
subrange. This problem will t»e discussed further in the next chapter. 
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Chapter 111 
SUBGRID SCALE M0DE13 


1 . 'n»e SCS Reynolda Streas 

In the preceding chapter we saw that there are terma in the equations of 
large eddy simulation that involve the amall or subgrid scale component o£ the 
velocity field, and, as this small-scale velocity field will not be computed, 
these terms must be modeled. This chapter will be devoted to a discussion of 
the models used for the so-called subgrid scale (SGS) terms. 

To begin, it is well to look at the physical significance of the SGS 
terms. Equations (2.2) and (2.3) describe the development of the large 
eddies. In them, the terms containing the small scale velocity represent the 
interactions between the large and small eddies. On the average, kinetic 
energy is transferred from the large eddies to the small ones, there is energy 
flow in both directions, but the net flow is usually toward the small scales, 
leslie and Quarinl (1979) estimated that the gross transfer to the small 
scales is about 1.5 times the net transfer. In other words, approximately 
one-third of the energy transferred to the small scales is returned. We shall 
see later that; the net energy flow may be in the reverse direction in some 
cases. The st bgrld scale terms in Eqs. (2.2) and (2.3) must represent the 
effect of these transfers on the large scales. In the normal situation, the 
net energy transfer to the small eddies appears to be a dissipation to the 
large eddies — energy lost that will not reappear. Thus the model should 
normally be dissipative. 

llie terras which need to be modeled were derived in the previous chapter 
and can be written; 


"ij - Vj-Vj 


“i“j “i“j ‘‘‘ “i“j 


(3.1) 


As we showed, we prefer to work with the SGS Reynolds stress defined by 


“ Si ' n \kSi 


(3.2) 
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It is ?Xbo worth mentioning at this point that the terms we have called 
the Leonard stresses; 


1 


ij 




(3.3) 


(which were first discussed by Leonard (1973)) may need special treatment. 
These terms are zero in the Deardorff-Schumann approach but not in the filter-' 
Ing method. Investigation has revealed that they are responsible for only a 
small amount of energy transfer between the large and small scales. Their 
major effect seems to be redistribution of energy among the various large 
scales . 

The contents of this chapter are as follows. In the next section, equa- 
tions governing the SGS Reynolds stresses will be derived and discussed. We 
shall also compare SGS modeling to Reynolds-averaged modeling. In Section 3, 
a computational method for validating SGS models will be described and some 
results given. This will be followed in Section 4 by a discussion of eddy 
viscosity models, the ones in most common use today. Section 5 will describe 
some of the contributions that theory has made to the state of the art in SGS 
modeling. Some new ideas about SGS modeling form the subject of Section 6. 
Higher-order modeling will be taken up in Section 7. Finally, Section 8 will 
discuss some effects that arise, when there are extra rates of strain (in Brad- 
shaw's sense) in the flow. We shall end the chapter with a short summary of 
the principal points. 


2. The SGS Stress Equations 

It is not difficult to derive a set of equations describing the dynamical 
behavior of the quantities R^^ defined by Eq. (3.1). However, the process 
is somewhat tedious. One takes the Navier-Stokes equations for u^ and also 
writes them with i replaced by j. The equation for u^ is multiplied by 
Uj and vice versa. Adding the two resulting equations and filtering the 


result yields an equation for 
the dynamical equation for u 


UfUj. 


By repeating the same procedure using 


one can derive an equation for 


Subtracting these two equations, we have the desired equation for R 


UfUj. 


ij* 
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t ^ “k ‘‘ij 
(Convection) 




9 u 


J 


+ (R4. + KJ 
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Ik' " '“jk " ''jk^ 
(Production) 




(Redistribution) 


- 2 


9u. &u 9u 3u. 

_i JL - _ — i — 1 

Sxj^S^ 

(Dissipation) 


( 3 ^ 4 ) 


+ Diffusion terms 


There are many diffusion termsi they are not written explicitly, as we ahali 
not need them. Here, is the Leonard stress defined by Cq. (3.3). An 
equatloi- for the subgrid scale turbulent kinetic energy, can be derived 
by Caking the trace of Cq. (3.4). Subtracting 6^^ times Che resulting 
equation from Cq . (3.4) gives an equation for 

All of the terms in Cq. (3.4) are analogous to terms in the familiar Key*- 
nolds stress equations of time-average modeling. The interpretations are also 
similar. However, the dlfferenci?s are quite important. Cqs. (3.4) contain 
more terms than the equations folc the time-average Reynolds stresses because 
somo items Chat are zero in time-average approach are not zero when filtering 
is used. In particular, note the appearance of the Leonard stress in Che 
production term and, more Importantly, the fact that Che production term is 
filtered. All of Che terms 1.;^ Cq. (3.4) can be computed by the methods de- 
scribed in the next section and models for them sC< 2 died, but this has not been 
done to date. 

The most common assumption in turbulence modeling is Chat production and 
dissipation terms dominate the turbulence budget, and, as a first approxima- 
tion, we can equate them and ignore the other terms. For the time-average 
equations, this approximation is reasonable when applied to the turbulent 
kinetic energy budget far from solid boundaries, but it is less v^lid for the 
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coaponent tiquationii because the redistribution terns nay be quite large. Near 
walls, the diffusion terms becoote quite iaportant and the approxiaation is 
even aore questionable. The low Reynolds nuabers in this region nay also 
affect the structure of the turbulence. Nonetheless, the "production equals 
dis8i>'ation argument" is frequently invoked. 

For LHS, the situation is soaewhat different. It is important to note is 
that the model is assumed to represent a local spatial average of the local 
instantaneous small-scale turbulence. This is quite different from what is 
modeled in time- or ensemble-average modeling and our understanding of subgrid 
scale turbulence (and consequently, our ability to model It) is more 
limited. This is compensated for by the fact that a large eddy simulation 
calculation of a given flow is less sensitive to modeling errors than is a 
Reynolds-averaged calculation of the same flow. 

In particular, because the small scales of turbulence are highly Inter- 
mittent, we expect gradients of subgrid scale quantities to be relatively 
large. If this is the case, it is probable that the convection and diffusioiy 
terras, which are ignored in many time-averaged models, are more important in 
SGS modeling. On the other hand, we have recently found evidence that the 
pressure fluctuations and, more particularly, the pressure-strain correlations 
reside mainly in the large scales (this will be presented in Chapter 3), and 
they may b«i lest important in SGS modeling than they are in conventional mod- 
eling. Despite these differences, most SGS models to date have relied on 
ideas developed for time average models. 

3. Computational Validation of SGS Models 

Two approaches are commonly used for developing and testing time-average 
models. One method, favored by Lumley, Reynolds and others, uses simple tur- 
bulent flows (usually homogeneous flows) to test the validity of the models 
and to determine the adjustable parameters. Tlie major objection to this 
approach is that the structure of homogeneous flows differs considerably from 
the flows one really wishes to simulate, and the constants may not be valid in 
more complex flows. Ihe other method, used by Spalding, Launder, and others, 
adjusts the parameters to fit flows similar to the ones that one wishes to 
calculate. This is difficult because many of the parameters must be adjusted 
simultaneously and this can be a difficult procedure. 
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It l8 even note dUIficult to develop ■odels lor the eubgrld ecalee. Ueta 
on the anall ecalee of turbulence ate quite icarce, and direct validation of a 
model ueing experimental data ie nearly inpoiaible. Confe<3qvient.ly , the con- 
atanta have to be found by other aetho/ia. One approach ia almoat completed 
baaed on theory and uaea the propertiea of the inertial aubrange> illly 
(1%7) and others have ahown that the constant in the model can be derived on 
this basia. Unfortunately, it is not aluaya poaaible to to assure that in a 
computation the cutoff between the large and small scales will lie in tlte sub- 
range, so one needs to be cautious about adopting the results of this 
approach. Indeed, a number of authors found it necessary to modify tlie SCS 
moiiol constant to obtain good results. 

There is a second approach. With the current generation of computers it 
is possible to compute homogeneous turbulent flows with no approximations 
other than those present in any numerical simulation. At present, it is 
possible to do such calculations with grids as large as 64 x 64 x 64 and, in a 
limited number of cases, 12B x 12B x 128. Ihis allows simulation at Reynolds 
numbers based on Taylor microscale up to approximately 40 (80 with the larger 
grids). The results can be regarded as realizations of physical flow fields 
and are an interesting and important complement to laboratory results. In 
particular, the computational results provide all three velocity components 
and the pressure at a large number of spatial points for a relatively short 
time span. Uie laboratory data typically give one or two velocity components 
over a longer time span at Just a few spatial points. 

liaving a realization of a flow, we proceed in much the same manner an 
experimentalist would. The computed field can be filtered to give its large 
scale component; the small scale component is obtained by difference. We can 
then compute the terms that need to be modeled, and, from the large scale 
field, we can also compute what the model predicts these terms to be. Direct 
comparison between the model and the exact value is then possible. This can 
be done in a couple of ways. 

One method is to use a scatter plot. The exact value of the SGS Reynolds 
stress at each mesh point is plotted against the value predicted by the model. 
If the model is correct, the results lie on a straight line; a totally invalid 
model produces a random pattern of points (usually a circle). This is a very 
graphic test of a model. Some scatter plots will be shown later. 
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Die eecond owthod le to coapare the aodel and exact reaulta ntatletlc- 
aliy. In ouc work we have need the correlation coefficient ae a aeaeure of 
the validity of a aodel. Diia ia a crude teat, but it aecao to be sufficient 
for our purposes. It is iaportant to recall that the square of tlie correla- 
tion coefficient ia approxlaately the fraction of the data that the aodel is 
correctly predicting. 

Diese are very severe tests of aodels — auch aore severe than the tests 
usually applied to Keynolds-averaged aodels. It is possible for a aodel which 
performs poorly in these teste to do mil in actual slaulatlons. tiowever, 
failure of a aodel to do well is a signal for caution. 

Use of this kind of testing for Reynolds-averaged aodels will be taken up 
in Chapter 5. 

4. Eddy Viscosity Models 

Eddy viscosity aodels can be '’derived’' from the "production equals dis- 
sipation" arguaent discussed earlier. This is done in a number of places and 
need not be repeated here. For subgrid scale turbulence, the eddy viscosity 
model amounts to assuming that the subgrid scale Reynolds stress is propor- 
tional to the strain in the large scale flow: 

/a'u 3u A 

The eddy viscosity Vj has the dimensions of a kinematic viscosity. 
Most work is leased on the at..ufflpCion that the eddy viscosity could be repre- 
sented by: 


- (CA)^ IS) (3.6) 

where A is the width associated with the filter and |U| ■ 

Recently, a number of authors have shown that this is correct only if the in- 
tegral scale of the turbulence is smaller than A . Since UES is designed for 
this not to be the case, it is better to assume that: 

- c A^^^ I’S:) (3.7) 
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where L Is the Integral scale o£ the turbulence. Usually L is estimated 
from L - where e is the dissipation. 

Eddy viscosity models have a long record of reasonable success In time- 
average modeling of simple shear flows, and one might expect them to do well 
as SGS models. In fact, they have been found to do well In some of the homo- 
geneout. flows^ In particular, for the homogeneous flows In which there Is no 
mean strain, one is able to predict most of the low-order statistical quanti- 
ties (for example, the mean square velocity fluctuations and spectra) quite 
well using eddy viscosity models; that the higher-order statistics, which are 
sensitive to the small scales, are not well predicted should be no surprise. 
In the homogeneous flows with strain or shear, there is evidence (McMillan et 
al. (1980), Shlrani et al. (1981)) that the energy transfer can be reversed 
and flow from the small scales to the large ones. In such cases, the model 
should no longer dissipate the energy of the large scales. Eddy viscosity 
models, which are guaranteed to dissipate energy from the large scales, cannot 
predict this behavior. Despite this, they may not function badly In actual 
slmulvHtlons. The reason Is that the smallest scales of the resolved field, 
from which the model normally extracts energy, become relatively weak In these 
flows, and the model may actually dissipate very little energy. Furthermore, 
the principal difficulty In computing these flows usually arises from the 
delivery of a significant amount of energy to scales larger than the compu- 
tational domain. Tills makes the normally used periodic boundary conditions 
incorrect, and the results cannot be relied upon. 

Eddy viscosity models are Incapable of handling other classes of flows. 
For example, In transitional flows, we must expect that most of the energy 
will be In the large scales i.e., the small scales are not In equilibrium with 
the large scales and the "production equals dissipation" argument Is 
Incorrect. Furthermore, although Moln et al. (1978) had reasonable success in 
simulating channel flow with these models, later extensions by ICim and Moln 
(1979) and Moln and Kim (1981) clearly show the deficiencies of the model. 
Th^ found that eddy viscosity models (several were tried) were unable to 
maintain the energy of the turbulence. Tlie problem Is only partially due to 
the model, as the turbulence tends to decay even when the model Is 
eliminated. This will be discussed in more detail In Chapter 7. 
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CXark et al. (1979) and McMillan and Kerzlgcc (1979), and Mcillllan at al. 
(19BU) have applied the model-teating method deacrlbed above to eddy viscoBity 
SOB niodclB. A typical scatter plot is shown In Fig. 3.1, in which the exact 
subgcld scale stress is plotted against the Smagorinsky model value. It can 
be seen that there Is a little correlation between the two data sets (the 
correlation coel ticlent Is approximately .4 for the case shown), but It is 
even more clear that this is far I'rom an adequate model. llxls result Is 
fairly typical, although there are variations in the correlation coefficent 
with many of the significant parameters. 

The results show that eddy viscosity models are rather poor and, in fact, 
they become even poorer wl»en there Is mean strain and/or shear In the flow. 
However, it is not easy to find more accurate models (we shall look at this 
below), so we may be forced to use eddy viscosity models until Homethlng 
better is developed. Furthermore, as McMillan and Ferziger have shown, the 
method can be used to predict the effect of Keynolds number on the model 
parameter. T1>clr results are shown in Fig. 3.2. When these results were 
applied to channel flow by Moin and Kim (private communication) , they did not 
produce the desired effects, probably for the reasons given above. 

In the above, we have used the fact that the natural length scale of the 
SGB eddies is the width, A , associated with the filter. By definition, this 
is the scale that defines whether an eddy is large or small, and there is 
little reason to suspect that this is not a correct choice. 

However, when the filter is anisotropic, as it should be in computing 
shear flows, it is not quite so clear what is the correct length scale. 
Almost everyone has used the cube root of the filter volume: 

A - 

However, Bardina et al. (1980) showed that a better choice might be: 

A - (Aj + d- A|J)W2 (3.9) 

It is recommended chat Eq. (3.9) be adopted for general use. 
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'Oie RoXe of 'Rieory 

Tlieorei..tcal Insight plays a considerable role In understanding the phys- 
ics of turbulence and contributes considerably to modeling Ito Tlirbulence Is, 
however, a problem of such complexity that the role of theory In our present 
etate of knowledge is smaller than In most areas of physics or engineering. 
Progress has been f rustratlngly slow. A review of recent theories Is given by 
Leslie (1973). 

Most theories provide llmltei Information about turbulence. Usually, the 
theories were developed for homogeneous turbulence and have proved difficult 
to generalize. 

The theories which have attracted the most attention are Rralchnan's 
direct interaction approximation and others related to it. lliese theories are 
statistical in nature i.e., they attempt to make statements about averages of 
turbulence quantities rather than the detailed dynamics, Ttie question of 
whether this theory could be extended so as to yield Information about the 
small scales of turbulence and, thus, to provide a SGS model has been inves- 
tigated by Leslie and his co-workers. 

Tlie theory necessarily deals with statistically averaged SGS turbulence. 
We imagine an ensemble of flows which have the same large-scale motions but 
different small-scale motions and ask for the average behavior of the small- 
scale motions. Whether this is adequate for modeling purposes is an open 
question, but the information generated should be helpful. 11118 theory, like 
many others, is capable of predicting the existence of an Inertial subrange, 
but, unlike most others, it can pred.vct the Kolomogorov constant as well. 

Love and Leslie (197b) extended the theory and showed that a form of the 
eddy viscosity model could be deduced from it. In particular, they predicted 
the constant in the model and showed that the large scale strain rate that 
appears in the eddy viscosity model ought not to be the local one but a spa- 
tial average. Hie constant predicted in this way Is in good agreement with 
that obtained by other theoretical arguments and from empirical fits to 
experimental data. 

With respect to spatial averaging of the strain rate in the eddy viscos- 
ity, the evidence is mixed. Love and Leslie (I97b) found that it was impor- 
tant in the solution of Burgers' equation, but Mansour et al. (1978) found 
that it did not matter much. 
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A number of other Issues were investigated by Leslie and Quarlnl (1979). 
In particular, they divided the SG3 terms into "outscatter” and "backscatter" 
terms representing, respectively, the energy flows to and from the subgrid 
scale. They found that eddy viscosity models appear to represent the outscat- 
ter fairly well, but they could not say much about the backscatter. 

Although limited, these theories are proving useful in choosing and vali- 
dating models. 


6. A Scale Similarity Model 

All models, by definition, relate the SGS Reynolds stress to the large 
scale flow field. Eddy viscosity models view as a stress and make an 

analogy between ir and the viscous stress. These models are guaranteed to ex- 
tract energy from the large scale field (i.e., they are dissipative). It is 
difficult to construct other models with this property. However, as noted 
above, the desirability of this property is questionable ic sheared and 
strained turbulence. 

It is important to observe that the interaction between the large and 
small scale components of the flow field takes place mainly between the seg- 
ment of each that is most like the other. The major interaction is thus be- 
tween the smallest scales of the large scale field and the largest scales of 
the small scale field (regions 1 and 2 of Fig. 3.3). This is what the SGS 
term in the filtered equations represents. Since the interacting components 
are very much alike, it seems natural to have the model reflect this. To do 
this requires that we find some way of defining the small scale component of 
the large scale field TT^ . One way to do this was suggested by Bardina et al. 
(1980). Since represents the large scale component of the field, fil- 

tering it again produces a field (u^^) whose content is still richer in the 
largest scales. Thus, 

\ “ “i ■ “j (3.10) 

is a field which cr/iitalns the smallest scales of the large scale component of 
the flow field. This suggests that a reasonable model might be 


T 

Ij 


CU^Uj 


(3.11) 
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or, better yet 


•ij 


(“i“j ■ “i“j) 


(3.12) 


This modification is sugijested by considering the “cross-terms," e.g., 

Preliminary tests have shown that this model is not dissipative, but it 
does correlate very well with the exact stress, a scatter plot is given in 
Pig. 3.4. Tills suggests that a combination of the two models might be better 
yet. Tlie correlation is largely due to the fact that, with a Uausslan filter, 
the two fields in question contain much the same structures. With other 
filters, particularly one which is a sharp cut-off in Fourier space, the 
correlation is smaller. These models are currently being investigated. 


7 . Higher-Order Models 

The inadequacies of algebraic eddy viscosity models in Reynolds-averaged 
modeling have been known for a long time. A number of mure complex models 
have been proposed, and, since they have analogs in SGS modeling, a brief 
review of them is in order. We shall go into some of these models in more 
detail later. 

Many of the improvements are based on the notion that proportionality 
between Reynolds stress and mean strain rate is valid, but the eddy viscosity 
formulation needs improvement. In these models one writes: 


- C^qil (3.13) 

where q and i are, respectively, velocity and length scales of the turbu- 
lence. In the simplest such models, the length scale is prescribed and a 
partial differential equation for the turbulence kinetic energy (q /2) is 
solved along with the equations for the mean flow field. These are called 
one-equation models; their record has not been particularly good, and most 
people now use still more complex models. In particular, the assumption of a 
prescribed length scale has been questioned, and methods of predicting the 
length scale have been proposed. Of these, the most widely used models are 
those in which an equation for the dissipation of turbulent kinetic energy 
(which really represents energy transferred to the small, dissipating eddies) 
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is added to the equations used in one-equation models. Hie length scale is 
related to the dissipation e by: 

C - C^q^/«, (3.14) 

and we have the so-called two-equation models. Hiis is the most popular 
method of computing time-average flow fields at present. 

Finally, the most recent development has been the use of Che full Rey- 
nolds stress equations. In two dimensions, three PUb's are needed to define 
the Reynolds stress, while in three dimensions, six are required. Clearly, 
this is a rather expensive approach. 

A way of avoiding Che computational cost of full Reynolds stress methods 
is obtained by noting that the convective and diffusive terms can frequently 
be neglected. If they are, and approximations are made to the redistribution 
terms, the equations reduce to algebraic ones. Algebraic models have become 
popular in recent years. However, there is doubt as to whether the neglect of 
diffusion is correct near the wall. 

All of these models have analogs in SGS modeling, and a number of them 
have been used. Let us consider them in the order in which they were intro- 
duced above. 

First, consider one- and two-equation models. They have as their funda- 
mental basis the proportionality of the SGS Reynolds stress and the large- 
scale stress. We saw earlier that Che Smagorinsky model (an algebraic eddy 
viscosity model) correlates poorly with the exact SGS Reynolds stress. Clark 
et al. (1979) looked at the behavior of one-equation models as well as an 
"optimized" eddy viscosity model. In Che latter, Che eddy viscosity was cho- 
sen, at every point in the flow, to give the best local correlation between 
the SGS Reynolds stress and Che large-scale strain. By definition, no eddy 
viscosity model can do better Chan this. It was found that Che correlation 
coefficient improved somewhat relative to the Smagorinsky model (from approxi- 
mately .35 to .50 in a typical case), but this still leaves the model far 
short of what we would like to have. The lack of correlation seems to be due 
to the difference between the principal axes of the two tensors. Hils model- 
ing assumption needs Co be changed if further improvement is to be obtained 
(cf. McMillan and Ferzlger (1980)), and more ccMnplex models are required. 
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Schumann (1973) also used one-equaCion modols without finding Improvement over 
algebraic eddy viscosity models. 

Next, recall the earlier remark that convection and diffusion are likely 
to be more Important In SGS modeling than they are in time-average modeling. 
This means that the approximations needed to reduce the full Reynolds stress 
equations to algebraic model equations are less likely to be valid In the SGS 
case. However, several authors have used algebraic models. The applications 
have been almost exclusively to meteorological and environmental flows in 
which stratification and buoyancy effects are Important. These flows are sen- 
sitive to small variations In both properties and model, making It difficult 
to assess the accuracy of a model with precision. To our knowledge, no 

applications of these models to engineering flows have yet been made. 

It is probable that, to obtain a significant Improvement over the Smag- 
orinsky eddy viscosity model, we shall need to go to full Paynolds stress 
models. This, of course, Is not something to be looked forward, to as the 
computing cost Is likely to be more than doubled. The only use of these 
equations to date was In meteorological flows by Deardorff (1972, 1973a, b), 
who reported a computer time Increase of a factor of 2.5. Furthermore, the 
results were not improved to the degree that he had hoped for. Although this 
Is discouraging, Deardorff 's simulation was considerably ahead of Its time and 
had the additional difficulties associated with buoyancy, so it is hard to 
make definitive conclusions. Thus, we cannot conclude much about these models 
at present, and quite a bit of work nesds to be done on them before they 
become useful tools of the trade. 

8. Other Physical Effects 

The author's group has done full simulations of the effects of compres- 
sibility on turbulence and the mixing of passive scalars In turbulent flows. 
To date, the work has concentrated on evaluating time-average models, because 
it was felt that this is the area In which the work has the most immediate 
impact . 

The effects of compressibility on SGS turbulence are probably quite 
small. The effect on the turbulence as a whole have been found to be fairly 
weak, except for effects due to the propagation of acoustic pressure waves. 
Since the latter are large scale phenomena and the Mach number of the SGS 
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turbulence Is small, we expect that compressibility will have only a weak, 
effect on SCS modeling. 

On the other hand, SGS modeling of turbulent mixing Is quite Important. 
If we are to simulate combusting flows. It will be necessary to treat the 
small scales accurately, since that Is where the action Is In these flows. 
Tlie effect of the Prandtl/ Schmidt number on time-average models Is modetrately 
strong, and we expect Its effect on SGS models to be even stronger. 
Furthermore, the specific effects due to combustion are also likely to be 
important on the small scales. The author Intends to look at SGS modeling of 
mixing and combusting flows In the near future. 

Another effect of considerable Importance in application is buoyancy, 
which was mentioned earlier in connection with the meteorological simulations. 
Flows in which buoyancy Is Important and, particularly, chose which are driven 
by buoyancy are very difficult flows Co measure or simulate, and a great deal 
of work will need to be done in this area. Important work in this area has 
been done by the Karlsruhe group (GroCzbach et al. (1979)), and further work 
is under way in London (Leslie (1980)). 

Finally, we should state Chat meteorologists and environmental engineers 
have a great interest in both mixing and buoyancy effects, and considerable 
effort in these areas has been made by these people. In particular, we note 
again the work of l>eardorff cited above and that of Sommeria (1976), Schemm 
and Lipps (1978), and Findikakis (1981). One of the principal difficulties of 
these flows Is Chat Che Reynolds numbers are so large that eddies of length 
scale equal to the grid size are quite Important. Consequently, Che SGS 
eddies do not behave entirely like "small eddies^” they carry a significant 
fraction of the total energy and are therefore hard to model. 

9. Summary of the State of SGS Modeling 

From the arguments given above, we can reach Che following conclusions 
about the current state of Che art In SGS modeling. 

1. Although they are Inadequate In detail, eddy viscosity models can be 
used In simulating homogeneous turbulent flows. However, they seem Co be 
inadequate for inhomogeneous flows, especially those In which solid boundaries 
are Important. 
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2. For models in which the length scalf; Is prescribed, the length scale 
of Eq. (3.9) Is preferred. 

3. One- and two-equation turbulence models are unlikely to provide 
significant Improvement relative to algebraic eddy viscosity models. An 
exception to this might be transitional flows. 

4. Full Reynolds stress models offer promise as future SGS models. 
However, the modeling assumptions probably need to be different from those 
used In time-average modeling. 

5. The scale-similarity model Is promising, but only when used In con- 
junction with other, dissipative, models. 

6. Full simulations seem to be the befit way available at present for 
testing SGS models and determining the parameters in them. Turbulence theo- 
ries can also be profitably used In this regard. 

7. Full simulations and large eddy simulations can both be used In time- 
averaged model building. Tills Is the area In which both types of simulations 
will make their greatest Impact in practical engineering calculations in the 
near future. 
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Chapter IV 
NUMERICAL METHODS 

1 . btathematlcal Preliminaries 

This chapter Is devoted to setting out the numerical methods used In full 
and large eddy simulations. To some extent, numerical methods are always 
tailored to the problem, higher-level simulations of turbulent flows are no 
exception. 

The partial differential equations governing a flow were given In Chapter 

2. To complete the methematlcal setting. It Is necessary to specify Initial 

and boundary conditions. This Is not easy. Higher- level simulations need 
details of the Initial state, and experimentalists are unable to provide suf- 
ficient data about the Initial state of their flows; some of the Initial data 
therefore has to be Invented. An equally serious problem is that, as the 
{^vler-Stokes equations are nonlinear. It Is not always known what boundary 
conditions should be specified, l.e., we may not know whether a problem Is 
well-posed or not. There are a number of examples of people attempting to 

solve mathematically Ill-posed problems. Another Issue Is that the partial 
differential equations have several conservation properties, and It Is Impor- 
tant that they be preserved In the numerical treatment of the problem. 
Finally, there are the difficulties Inherent In the numerical methods them- 
selves — accuracy, stability, and aliasing, among others. All of these, need to 
be considered. 

The equations governing incompressible flows are of mixed type; they 
contain elements of both parabolic and elliptic parclal differential 
equations. This is a consequence of the momentum equations containing time 
derivatives, but the continuity equation not having any. As a result, one 
cannot advance the continuity equation In time. These equations are called 
Incompletely parabolic by mathematicians. Means of dealing with both types of 
behavior are needed. The compressible equations, which are hyperbolic, are 
actually easier to deal with from a numerical point of view. 

All of these issues will be taken up in the remainder of this chapter. 
Additionally, we shall need to describe the numerical approximatlcvis used in 
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the computations. Throughout the ctiapter, It is Important to keep In mind the 
kinds o£ flows that we are trying to simulate. Tliey are geometrically simple 
turbulent flows. Die fact that they are turbulent means that the high wave- 
number components of the velocity field are large. Large gradients of the 
variables can occur in any part of the flow, this has an important influence 
on the choice of numerical methods, Un the ocher hand, the simplicity of Che 
geometry helps considerably in developing accurate numerical methods. 

2. Boundary Conditions 

Die simplest flows to be simulated are the homogeneous turbulent flows. 
By definition, these flows are statistically identical at every point in the 
flow. For these flows, the most convenient and most accurate boundary condi- 
tions are perxodic ones. Die portion of the flow within a rectangular paral- 
lelepiped is simulated, and the boundary conditions prescribe that the state 
of the fluid at a point adjacent to any of the boundaries is identically that 
on the opposite face of the parallelepiped. These conditions avoid Che need 
for specifying the details of a highly chaotic motion on the surfaces and are 
the most realistic means of enforcing the idea that any point in the flow is 
indistinguishable from any other point. 

Diere is one point that requires extra care. In homogeneous turbulent 
flows on which mean straining or shearing flow is imposed, it is convenient to 
solve for Just the part of the flow containing the turbulent fluctuations; the 
mean flow is eliminated. When this is done, it is found that there are terms 
in the equations that uo not admit Che use of periodic boundary conditions. 
It is Chen necessary Co do the computation in a coordinate system that deforms 
with Che mean flow, and the ability to use periodic boundary conditions is 
restored. Diis will be taken up again in Chapter 5. 

The only other flows chat we shall consider in any detail in this report 
are inhomogeneous in one coordinate direction. Of course, this means that 
they are homogeneous in Che other two directions, and these directions can be 
treated by the periodic conditions described above. There are two types of 
conditions we must deal with in the inhomogeneous direction; they follow from 
the nature of the flows we shall be simulating. 
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For free shear flows, we would like to prescribe the condition that the 
flow la at rest infinitely far froo the shear layer. Dealing with an infinite 
region is difficult, and two methods have been used for this problem: 

i) One can use a finite computational domain. At the top and bottcMS of 
the domain one specifies that the vertical derivatives of the horizontal 
components of the velocity are zero, and the component of the velocity normal 
to the boundary is set to zero. These are known as no-stress boundary condi- 
tions. Unfortunately, no-stress conditions imply the existence of image flows 
outside the computational domain: the images are reflections of the flow in 
the boundaries. To a»>aure that the image flows do not interfere with the 
physical one, there must be no vorticlty close to the no-stress boundary. 
Tliis means that a considerable portion of the computational domain must be 
wasted in computing the potential part of the flow. 

li) One can use a coordinate transformation that maps the infinite do- 
main onto a finite one. Standard numerical methods can then be used. It is 
Important to choose a mapping that is compatible with the method used foe 
evaluating derivatives. Tills issue will be dealt with in more detail later. 

Ttie second type of inhomogeneous flow that we shall consider is fully 
developed turbulent channel flow. TVo different approaches have been taken 
for simulating this flow: 

i) Deardor ff and Schumann decided not to treat the wall directly. The 
reasons will be stated In detail In Chapter 7. Instead, they decided to com- 
pute only the part of the flow within and beyond the region In which the 
velocity profile is logarithmic. The boundary conditions must then assure 
that the velocity profile be logarithmic at the edge of the computational 
domain. In addition, it is necessary to specify something about the nature of 
the turbulent fluctuations at this boundary. They assume a relationship 
between the velocity and stress fluctuations at the boundary, this is the 
simplest assumption one can make, and there is no evidence for any other 
choice, but it has been called Into question. 

11) One can compute the entire flow. Including the region near the wall. 
The wall conditions are then the no-slip conditions that must be imposed at 
any solid surface. This is a much simpler boundary condition to deal with 
numerically. 'fhe price one pays is thiit all of the small structures near the 
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wall Bust be coaputed explicitly; this leads to considerable diificulty, as we 
shall see in Chapter 7. 

3. Treatment of the Spatial Derivatives.* Conservation Properties 

In all flow computations, the spatial derivatives are approximated in 
terms of the values of the dependent variables at grid points. Higher-level 
turbulence computations are no different from others in this respect, the 
methods used in these flows are also used in other types of flow simulation. 
Again we note that the geometric simplicity of the flows treated by higher- 
level simulations allows use of methods that might not be easily applied in 
more complex geometry. 

Before giving the specific approximations to be used, it is important to 
discuss conservation properties. We believe that this issue is not emphasised 
sufficiently in the literature. Tlie dynamical equations are essentially 
microscopic conservation equations. The continuity equation expresses 
conservation of mass. In the compressible case, the Havier-Stokes equations 
express momentum conservation (or what is the same thing, Newton's second 
law) , and there is a separate energy equation to express the fact that total 
»'nergy is conserved. In the incompressible case, the Navier~Stokes equations 
still express momentum conservation, but, in the absence of an explicit energy 
equation, they are also responsible for conserving the only significant energy 
in the flow--the kinetic energy. This leads to one of the principal 
difficulties in the treatment of incompressible flows. 

By integration of the microscopic conservation equations over a finite 
volume, we obtain macroscopic conservation equations. For the incompressible 
form of the continuity equation we obtain the global conservation of mass 
relation: 
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The Navler-Stokes equations give rise to the well-known momentum theorem: 
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Finally, multiplying the Navier-Stokea equations by u^ and integrating over 
a iinite volume, we obtain ttie equation of kinetic energy conservation: 
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Each of these equations states that the conserved property changes only by 
flow of the property through the bounding surface, this is a consequence of 
the fact thof there are no sources of any of these properties within the vol- 
ume. If periodic boundary conditions are applied, the surface terms integrate 
to zero. In Eqs. (4.1), (4.2) and (4.3), S is the surface of the volume V. 

The kinetic energy conservation equation (4.3) is especially interesting. 
The only non-surface term is the viscous dissipation term, which is usually 
small. It is essential to note that the kinetic energy within the control 
volume is not changed by the convection and pressure gradient terms and that 
the chain rule sov)’ ■ u'v + v'u) and the continuity equation are used in 
eliminating the volume integral of the pressure. 


It is crucial that the numerical approximations to the equations retain 
these properties. For the continuity and momentum conservation, this is usu- 
ally not difficult. It usually turns out that, if the equations are written 
in the proper form (the so-called conservative form we have used throughout), 
then almost any approximation will yield these conservation properties. The 
principal difficulty is with the kinetic energy. Normally, the verification 
that the numerical approximation guarantees energy conservation has to be done 
on a case-by-case basis. A means of avoiding this difficulty was found by 
Mansour et al. (1978). If the Navler-Stokes equations are written in the 
form: 
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rather than (2.1), the derivation of the conservation of kinetic energy equa- 
tion (4.3) can be based on a symmetry property, and the use of the chain rule 
can be avoided. Since numerical approximations do not always have a chain 
rule but the symmetry property always holds, using the Navier-Stokes equations 
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In the form (4.4) can simplify the job of finding well-behaved numerical 
methods. 


Many Workers (Deardorff (1970), Schumann (1973), Antonopoulos (1981), 

Shaanan et al. (1975), among others) have used the stagge red-grid mesh 

system. The grid Is shown In Fig. 2.2 for the two-dimensional case, the 

variables are given at the mesh points shown in the figure. The control 

volumes for the various equations are different and are displayed In the 

figure, We shall not give the finite difference equations here, as they appear 

In several other works. fhls grid system has the nice property that all of 

the conservation properties are obtained without difficulty, and, as we shall 

see In the next section. It gives no problem with the calculation of the 

pressure. It Is the natural grid system for the Incompressible equations and 

has been used more widely than any other. Part of the reason for the success 

of the stagger ed me sh system was explained by Shaanan et al. (1975). The 

approximation u^u . u^u. which has been used by Deardorff and Schumann Is 
1 j 1 j 

valid In the staggered grid system, because the truncation errors represent 
the difference between c;'.ese two terms (the Le onard stress) quite well. 

Stated otherwise, the staggered grid approximates <nore accurately than 

It does thus leads to great simplification In the finite difference 

equations . 


If a re<^^r;lar grid Is used. It Is necessary to use a fourth-order finite 
difference method In order to assure that the i.eonard stress Is properly com- 
puted. This can be done, but the method is cumbersome (Kwak et al. (1975)). 


Another popular method of computing derivatives In directions In which a 
flow Is homogeneous Is by means of Fourier transforms — the pseudospectral 
method. In this method one uses the discrete Fourier transform. Any function 
defined at a set of equally spaced mesh points x^ ■ jAx, j ■ 1,2,...,N, can 
be represented by the discrete Fourier series: 
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which differs from li)q. (4.5) only by the sign of the exponent and the factor 
i/N; thuSi bpoth transforms can be computed In the same way. these results 
can be used In the following way: Given the values of the function f(x) on 

A 

the grid points Xj ■ JAx, we can compute f(k^£) from hq. (4.6). When theue 
are used in Eq. (4.5) and Xj Is replaced by the continuous variable Xt the 
result Is an interpolation formula. As such, It can be differentiated with 
respect to x, and this provides a method of computing spatial derivatives. 
In fact, specializing the result to the grid points, we have: 
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The derivative df/dx can be computed by using the discrete values f(xi) to 
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compute f(k^), multiplying the result by ^od computing the Inverse 

transform. 'Hie result is an extremely accurate estimate of the derivative. 
This method is especially well adapted to the calculation of the derivatives 
of periodic functions, which explains its widespread application in the compu- 
tation of homogeneous turbulent flows. 


llie practical use of the Fourier trans.form as a numerical tool is made 
possible by the existence of an extremely fact algoritlim for its computation — 
the so-called fast Fourier transform (FFT) algorithm. 


For Later application it is important to note that this method could also 
be used to compute finite differences. It is not effective to use this as a 
tool for computing derivatives, but it can play an important role when we come 
to solving the equation for the pressure. As an example, we take the standard 
second-order central difference approximation: 
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Tlie derivative obtained from this formula can be put into the form of Eq , 


(4.7) with ikj 


replaced by Ik^ “ i(sin k^Ax)/Ax; we call kj^ 


the effec- 


tive wavenumber. Effective wavenumbers are a good way to measure the accuracy 
of finite difference methods that are required to differentiate functions 
which contain significant high-wavenUvmber components, and it is not difficult 
to derive the effective wavenumber for various finite difference approxima- 
tions. Some effective wavenumbers are plotted in Fig. 4.1. 
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Next, let U8 dlscuBs the treatment of directions in which the flow is not 
homogeneous. For the free shear layer, we noted in the previous section that 
there are two ways of dealing with the direction normal to the flow (the shear 
direction). When the no-stress boundary conditions are used, one can general- 
ize Che Fourier method described above. The key idea is to expand the func- 
tions in terms of sines or cosines (using the set appropriate to the boundary 
conditions for the particular function) rather Chan the complex exponentials 
of Fq. (4.5). The I'lumerical algoritlim for computing sine and cosine trans- 
forms is equivalent Co computing the exponential transform (4.5) using 2N 
rather than N points. Tljus the cost of computliig the derivatives is approx- 
imately doubled when this method is used. We noted earlier that this approach 
suffers from loss of accuracy due to images. 

The alternative to the use of no-stress conditions is the use of a trans- 
formation which Cakes Che physical coordinate z into a computational coordi- 
nate C : 


z " h(c) 


(4.9) 


such that - ® < z < “> transforms to - 1 < 5 < 1. The derivative becomes: 


. i_ -1- 

dz h' dC 


(4.10) 


Tlie trick to making this a successful method is to choose the transformation 
such that l/h' can be expressed in terms of just a few low-order sines 

and/or cosines. It is Chen possible to obtain accuracy almost as good as that 
of Che Fourier method for infinite regions. Details of this method are given 
in the report by Cain et al. (1981). 

Finally we come to dealing with directions in which there are solid 
walls, i.e., a numerical method for Creating channel flow. Tliere are two 
choices that have been commonly used. The first is the use of Chebychev poly- 
nomial expansions. Tills is equivalent to a Fourier method on a nonuniform 
grid and has been used by Kells and Orszag (1980) and by Orszag and Patera 
(1980); see also Kim and Moin (1980). 

The other method for creating the channel is to use a finite difference 
method on a nonunlfonn grid; this is equivalent to using a coordinate trans- 
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formation in this direction. llie choice lias generally been adopted (Moin at 
al. (1978), Moin and Kim (1981)). We shall deal with this further in Chapter 

7. 

Another important Iscue is aliasing. Aliasing is the error introduced 
when two Fourier waves are multlpllodi this happens implicitly when the non- 
linear convective terms are computed. The waves resulting from the product of 
two Fourier waves contains the sum and/or difference of the original wavenum- 
bers, Tliese may fall outside the range of wavenumbers (-a/A < k < n/A) 
which can be carried in the calculation. When this liappens, the wavenumber 
which falls outside computational range is misinterpreted ("aliased") as one 
of the wavenumbers which lies inside the band. The result is a numerical 
error which, in mild cases, adds to the normal truncation error of a finite 
difference approximation and, in severe cases, can cause the calculation to 
become totally Inaccurate or even unstable. 

Aliasing can be controlled in two ways. The simplest way is to assure 
that the high wavenumbers are relatively unpopulated. Since these are the 
ones that cause the problem, eliminating them also solves the problem. In 
large eddy simulation, one can assure that the high wavenumbers are relatively 
unpopulated by using a filter which cuts off at a moderate wavenumber. In 
full simulations, the best way to control the problem is to keep the Reynolds 
number low. 

Tl\e other method of controlling aliasing is to compute the portion of the 
field which will be aliased and explicitly eliminate :! ;. This requires extra 
computation, but It allows one to include more energy in the high wavenumbers 
and the extra resolution gained may be worth the cost. 

4 . Time Advancement 

We now come to the method of advancing the solution in time» One of the 
first issues that arises is that of selecting an explicit method or an impli- 
cit one. It is important to remember that in higher-level simulations one is 
looking for time-accurate solutions to the equations of motion. This con- 
trasts strongly with relaxation methods, in which the object is to reach 
steady state as quickly as possible. The point of view that we adopt is that 
a well-balanced, time-accurate method is one in which the errors caused by the 
time advancement method approximately equal those introduced by the spatial 
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Introduced by the spatial differencing method. Once spatial difference ap- 
proximations and the time-advancement method have been chosen, this criterion 
selects the time step. The time advancement method must be stable for the 
time step' so chosen. It Is usually the case that the time step found In this 
way Is well within the stability bounds of explicit methods, so there Is no 
need to pay the extra cost associated with an Implicit method. Thus, with a 
few exceptions, noted later, the time-advancement methods used In higher-level 
simulations are explicit. -The common choices have been second-order methods 
such as leapfrog and Adams-Bashforth and the fourth-order Runge-Kutta method. 
Tiiese are standard methods of numerical analysis, so the formulas will not be 
given here. 


For purposes of discussing time-advancement methods. It Is convenient to 
rewrite the Navler-Stokes equations In the form: 
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dt 
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where the viscous and convective terms have been Included In Tliere Is no 
difficulty In time-advancing this equation by an explicit method. Most of the 
difficulties in solving the Incompressible equations come from the lack of a 
time derivative in the continuity equation; the compressible equations have no 
such problem. One method of avoiding this difficulty Is to treat the flow as 
if it were compressible and Iteratively drive the compressibility effects to 
zero; the Iterative nature of this process mak^';8 It Inefficient, however. 


A more efficient procedure is to note that application of the divergence 
operator and use of the continuity equation on Eq. (4.9) gives the Poisson 
equation for the preesure: 
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(4.10) 


When one looks at the time-discretized version of Eq. (4.9), it is found that 
forcing the pressure to satisfy the Poisson equation (4.10) at time step n 
guarantees that continuity will be maintained at time step n + 1. The mixed 
nature of the equations is brought into clear focus. The Navler-Stokes equa- 
tions (4.9) are treated as parabolic partial differential equations, but the 
pressure must be calculated from the Poisson equation (4.10), which Is ellip- 
tic. 
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One further Important point needs to be made here. Recall that, If the 
Navler-Stokes equations In the form (2.1) are used, then the derivation of the 
energy-conservation equation (4.3) requires use of the chain rule and the con- 
tinuity equation. If we are to have numerical energy conservation. It Is 
necessary to derive the numerical equivalent of Eq. (4.3). Assuming that ths 
required analog to the chain rule exists, the choice made for the numerical 
approximation to the pressure gradient dictates the numerical approximations 
used In the continuity equation. Otherwise, one cannot obtain energy conser- 
vation;, the usual consequence Is an unstabl'^ calculation. For example. If the 
central difference approximation Is used to estimate dp/9x. It must be used 
for the continuity equation as well. If a backward difference Is used for the 
pressure gradient, the continuity equation must use the forward difference 
oper'ator, and vice versa; this Is what is done on the staggered grid. 

Furthermore, one Is not free to finite difference the Poisson equation 
(4.10) arbitrarily. The correct approximation Is derived by applying the 
numerical divergence operator obtained in the manner described In the pre- 
ceding paragraph to the finite difference version of the Navler-Stokes equa- 
tions. Thus, the choice of the finite difference approximation for the 
pressure gradient dictates the method of differencing the Poisson equation. 
For example, if the central difference operator Is used for 3p/8x, it turns 
out that the difference operator for the Poisson equation must be the second- 
order central difference operator (as one might expect), but the grid spacing 
must be 2Ax and not Ax. We reiterate that the function of the Poisson 
equation Is to maintain continuity In the numerical sense; It is more Impor- 
tant to solve the correct equation than to obtain the most accurate solution 
to the exact partial differential equation. 

The most efficient method of solving the Poisson equation is by means of 
the fast Fourier transform. This Is the case whether one uses finite differ- 
ences or the pseudospectral method, the spatial derivatives. When finite 
differences are used, one can solve the Poisson equation by using Fourier 
transforms, but one must be careful to use the effective wavenumber rather 
than the exact wavenumber. 

The staggered grid method accomplishes all of this very efficiently. It 
does this so well that the need for being careful with finite difference meth- 
ods Is often overlooked. 
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There is one case In which we cannot use explicit methods. In the compu- 
tation of flows with solid boundaries, it is necessary to use a very find grid 
in the direction normal to the wall, close to the wall. A consequence is tliat 
the time step allowed by stability is then smaller than the time step allowed 
by the accuracy criterion. The principal difficulty comes from Che viscous 
term. In this case, it is necessary to treat Che viscous terms containing 
derivatives in Che normal direction implicitly. In fact, a special numerical 
method had to be Invented for this problem^ it will be sketched in Chapter 7 . 

5. Initial Conditions 

Tlie initial conditions for higher-level simulations cannot be derived 
directly from experimental results. The data never contain enough information 
to construct a complete initial field. In fact, the reported results of some 
experiments are quite incomplete and leave Che computor so much freedom that 
it is always possible to find initial conditions Chat allow the simulation Co 
match the experiment. From the point of view of one doing higher-level simu- 
lations, an ideal experiment reports not only the mean velocity and turbulence 
intensities, but information about the length scales as well. Ideally, com- 
plete spectral Information should be provided. 

We begin by considering the construction of a velocity field for the 
simulation of homogeneous isotropic turbulence^ the velocity fields required 
by the other cases are frequently derived from this. The task is to create 
an initial field that has a specified, energy spectrum and is divergence- 
freel There are several ways to do this; of these, Che following is one of 
the easiest. There are three steps in the process: 

1. Each component of the velocity at every grid point is assigned a 
random value. Tire resulting field is not divergence-free, nor does it have 
the desired spectrum. 

2. The curl of the field is taken; the resulting field is divergence- 
free. The numerical operator used to take the curl must be the same as the 
operator used to define the divergence. 

3. The Fourier transform of the velocity field is taken and in each 
Fourier mode is assigned an amplitude required to give the desired spectrum. 
The Fourier transform is Inverted, and the result is the desired initial 
field. 
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This procedure Is easily modified to give an initial field which is 
anisotropic. This can be done by biasing the random numbers used in the first 
step of the process. 

For flows in which there is a mean velocity profile (specifically, the 
mixing layer and the channel), it is necessary to give the mean velocity 
profile in addition to the turbulence. The method of producing an initial 
turbulence field must also be modified. In the case of the mixing layer, we 
want the fluctuations to be more Intense near the central plane of the flow 
than near the edges. Such a field can be produced in a manner similar to that 
described above. Insteady of allowing the field created in Step 1 to be uni- 
formly distributed in space, we give it the desired spatial distribution. The 
steps for removing the divergence and producing the given spec' rum are then 
essentially as described above. 

For the channel flow it was found that the subgrid scale model destroyed 
too much energ; and tended to make the flow become laminar if conditions of 
the kind described above were used. To prevent this, it was necessary to 
introduce large structures into the flow. These were obtained from solutions 
of the Orr-Sommerfeld equations. Although these are not the correct large 
structures for a fully developed channel flow, they are apparently similar 
enough to them. Randomness was added to the flow by introducing a small 
amount of more or less Isotropic turbulence which is divergence-free and Is 
zero at the walls. 
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Chapter V 

HOMOGENEOUS TURBULENCE 


^ • Claaslf Icatlon 

A homogeneous turbulent flow is one in which each point in the flow is^, 
in the statistical sense, equivalent to every other point. Ideally, this re*- 
quires an infinite medium of fluid, every part of which experiences the same 
forces. In practice, close approximations to these flows are produced in wind 
tunnels. The mean flow is designed into the tunnel, while the turbulence is 
usually created by a grid (or, in a few cases, by a set of Jets) and carefully 
controlled. The time evolution of the flow is simulated by observing its de- 
velopment as it moves downstre^un in the tunnel and invoking Taylor's hypothe- 
sis. If the gradients of mean quantities and other parameters of the flow are 
carefully chosen, an accurate approximation to a homogeneous flow is produced. 
It is not difficult to show that homogeneity requires the mean flow to be one 
in which the mean velocity is a linear function of all of the spatial coordi- 
nates. This severely limits the possibilities. 

Nearly all turbulent flows of engineering interest are inhomogeneous, the 
inhomogeneity is usually the result of the shear varying through the flow. 
When the Reynolds stresses in these flows are modeled, five separate effects 
are commonly considered. They were mentioned in Chapter 3 and are repeated 
here : 

a. Production . The creation of new Reynolds stresses via the inter- 

action of the Reynolds stresses with the mean flow. 

b. Dissipation . The destruction of turbulent energy and Reynolds 
stresses by the action of viscosity. 

c. Redistribution . The conversion of one component of the Reynolds 
stress into another without change of the total turbulent energy. 
Much of this effect is mediated by the pressure. 

d. Convection . The convection terms usually require no modeling, but 

their inclusion makes the local Reynolds stresses depend on the mean 
field in other parts of the flow. 

e. Diffusion . The carrying or Reynolds stress from one part of the flow 
to another via the self-interactions of the turbulence. 
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By definition^ homogeneous flows have no convection or diffusion) so we need 
to deal with, at most, production, dissipation, and redlstrlhutlon. 

Homogeneous turbulent flows can be grouped Into three categories accord- 
ing to the phenomena contained In them. Tlie first group contains the one flow 
in which the only Interesting effect is dissipation# (Inertial energy trans- 
fer among the wavenumber components is, of course, an element In all flows but 
Is not counted separately#) 

• Homogeneous Isotropic 'Rirbulence # This flow, which at one time was 
heavily studied because It was thought that It might provide the Insight 
Into the nature of all turbulent flows. Is the decay back to rest of 
fluid which has been set Into random motion. It Is still used as a means 
of finding turbulence model constants associated with dissipation and Is 
usually the first flow simulated by people doing higher-level 
simulations . 

• The second group of flows contains those In which there is exchange 
between the various components of the Reynolds stress (redistribution) In 
addition to dissipation, but there Is no direct production of turbulence 
energy. There are two such flows. 

a. Homogeneous TUrbulence with Rotation . The effect of rotation on 
isotropic turbulence is to produce anisotropy. Tlie effect Is primarily on the 
length scale and reduces the rate of decay of the turbulence. 

b. Return to Isotropy . Turbulence which has been made, anisotropic by 
the action of strain (see below) tends to return toward isotropy if the 
additional force is removed. 

• The final group contains the flows in which all of the phenomena that are 
possible in homogeneous flows actually occur. There are two major flows 
of this type. 

a. Strained Homogeneous Thrbulence . Turbulence which is Initially 
Isotropic (or nearly so) is put through a wind-tunnel section in which a fluid 
element is stretched in one direction and compressed in either one or two di- 
rections. The result is irrotatlonal strain which interacts with the existing 
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turbulence; there Is considerable turbulence energy production, and the flow 
becomes quite anisotropic. 

b. Sheared Homogeneous TVirbulence . Nearly Isotropic turbulence is 
produced in a flow which has uniform shear (a straight-line velocity profile). 
The effects are similar to those observed in the strained turbulence case. 

Tlie experimental data for these flows have been revi.ewe<‘i in a paper by 
the author (Ferziger (19BU)). 

We shall also consider flows with compressibility and mixing of a passive 
scalar. 

All of the flows described in this chapter arc one^/ which develop in 
time. It is uncertain that any of them reaches a steady state or even a self- 
similar state. Ihis issue is controversial; some authors believe that a self- 
similar state will be reached, while others do not believe so. In any case, 
these flows are sensitive to the initial conditions. In turn, this means that 
caution is required in interpreting them and that careful documentation of the 
initial conditions is necessary. 

All of these flows have been calculated by both full and large eddy 
simulation. Tlie results show that all of the physical phenomena observed in 
the laboratory have been shown to be a valuable tool in evaluating turbulence 
models. Much of the work in the area of model validation is recent and unpub- 
lished, and we shall give a brief overview of some of the principal results. 
It is also worth pointing out that a complete compendium of results from full 
simulations of homogeneous turbulence is being assembled by Dr. R. S. Kogallo 
of IilASA-Ames Research Center and will probably be available in the summer of 
1981. His results should be an important resource for people developing tur- 
bulence models. 

2. Isotropic Tlirbulence 

As we have mentioned earlier, isotropic turbulence is the simplest tur- 
bulent flow. It is therefore an obvious first target for any method of simu- 
lating or modeling turbulent flows. It has long been used by the developers 
of Reynolds-averaged models as a basis for choosing the con8tant(s) associated 
with the dissipation. It has also been a popular choice as the first flow to 
be simulated by higher-level methods, and it has been used extensively as a 
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basis ior testing subgrid'scale fflodels. We shal.’ review this work brielly in 
this section. 

To simulate these flows numerically, one begins with an initial condition 
that has the desired energy spectrum and is divergence-free. ^iethods of 
constructing such fields were described in the preceding chapter. In full 
simulations it is not necessary to begin the calculation with a realistic 
spectrum; one will develop in time. Of course, if one is trying to match an 
experiment, the experimental spectrum ought to be specified. In large eddy 
simulations of this flow, the initial spectrum is obtained by filtering the 
experimental spectrum. 

Ttie initial condition defines the initial Reynolds number. The Reynolds 
number commonly used to characterize this flow is based on the Taylor 
microscale X and the turbulence intensity q. Although these may not be the 
optimum choices, we shall follow custom and use them. In this flow the 
turbulence Intensity decays and the microscale increases with time, but the 
microscale Reynolds number decreases. 

At the first few time steps, the flow field cannot be regarded as repre- 
senting true turbulence. Tlie initial field does not contain the proper 
higher-order statistics or correlations; only after at least some of these 
have developed can the field he taken as representing physical reality. We 
have generally taken the behavior of the skewness of the velocity derivative: 

S - < (3u/3x)^ > / < Ou/8x)^ (5.1) 

as the measure of thi; quality of the flow field. It is nearly zero in the 
initial field and quickly rises to an asymptotic value at which it tends to 
remain for a considerable time, except at low Reynolds numbers. The time 
period in which the skewness Is rising is considered a "development” period. 
This is followed by a period in which the flow is realistically simulated. 
Finally, the size of the large structures grows to an appreciable fraction of 
the size of the computational domain, and periodic boundary conditions are no 
longer valid. At this point, the flow is no longer realistic, unphysical 
behavior is observed in the results, and the program has to be stopped. 

First, let us consider full simulations of this flow. Using 64 x 64 x 
64 mesh points in a calculation, one is able to compute at Reynolds numbers 
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up to - 50. This Is the practical limit on moat present computers, a lew 
cases have been run on a 12b x 12B x I2B grid which allows ttte Reynolds 
number doubled. Ihese Reynolds numbers are on the low end of the experimental 
ones; most experiments have been run with R^)^ in the range SO-AQO. The 
results of these computations match the experiment very well in terms of the 
decay of the turbulence intensity, the growth of the length scales, and the 
value of the skewness. Typical results are shown in Figs. 5. 1-3.3. 

The principal use to which these results have been put has been in the 
development and testing of subgrid scale models. Clark et al. (197V) and 
McMillan and Ferziger (1979,1980) used flow fields generated by full simula- 
tion of Isotropic turbulence in the way suggested in Section 3.3. Some of the 
principal results of this work were: that the Smagorinsky model correlates 
very poorly with the actual SGS Reynolds stress (the actual correlation coef- 
ficient is typically .30-. 40), that the width of the filter used in large eddy 
simulation ought to be at least twice the grid size, that changing the shape 
of the filter matters little, and that the model "constant" (which really 
ought to be called a parameter) is a function of Reynolds number that can be 
derived from this type of calculation. Since these results were covered in 
Chapter 3, we shall not repeat them here. 

Full simulation has also been used to study isotropic turbulence at low 
Reynolds number, a purpose for which it is ideally suited. At low Reynolds 
numbers (R^ < 10), it is possible to do full simulations with only lb x ib 
X 16 mesh points. Interest in these flows centers on the decay rate and the 
skewness. The decay of isotropic turbulence can be represented by: 

q2 - A(t-t^)"" (5.2) 

Theory shows that the decay exponent (n) is 2.5 at very low Reynolds number, 
and both theory and experiment show it to be approximately 1.2 at high Rey- 
nolds number. It is therefore of interest to compute the decay exponent as a 
function of Reynolds number. The results are compared with experiment in Fig. 

5.4. 

The velocity derivative skewness defined by Eq. (5.1) is approximately 
.5 at high Reynolds number and can be shown to drop to zero as the Reynolds 
number goes to zero. The direct simulation and experimental results are shown 
in Fig. 5.5. Figures 5.4 and 5.5 are from a report by Shirani et al. (19B1). 
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Now let U8 turn our attention to large eddy aifflulations ot this t'low* 
The major advantage of large eddy almulatlon la that, alnce the email eddies 
are modeled, the computation time is considerably reduced for a given Reynolds 
nuiuber. Alternatively, It is possible to go to higher Reynolds number with 
LES than with direct simulation. 

When large eddy simulations of homogeneous isotropic turbulence were 
first made, the results of full simulations were not available. Consequently, 
the constant had to be chosen to fit the decay of the turbulent kinetic 
energy. It was found that the same constant can be used whether Ib^ or 32^ 
mesh points were used; it was later found that the value obtained in this way 
agreed with those obtained from direct simulation to within 1U%. It is also 
in good agreement with theoretical estimates (Lilly (1967)) despite the fact 
that these flows are at Reynolds numbers too small to support an inertial 
subrange. Since the constant needs to be adjusted by this amount to account 
for changes in numerical method (mainly changes in the spatial differencing 
method) , this was one of the most important early successes or large eddy 
simulation. 

It was found that it made little difference whether the primitive Navier- 
Stokes equations or the vorticity form of those equal:) ons were used; it made 
very little difference whether the model was based on the strain rate or the 
vorticity; and it made very little difference which filter was used. However, 
if pseudospectral differencing is applied to the original Smagorinsky model, 
the shape of the spectrum at high wavenumbers is distorted. To remedy this 
problem, it was found necessary to evaluate all of the derivatives that occur 
in the model by second-order finite-difference approximations. This is simi- 
lar to the finding by Love and Leslie (1976) that the model ought to be aver- 
aged over a finite volume. A typical result obtained by large eddy simulation 
is shown in Fig. 5.6; other curves are similar and therefore not shown. 

Finally, it was found that large eddy simulation is incapable of comput- 
ing the higher-order statistical quantities such as the skewness and flatness 
with sufficient accuracy. These quantities are strongly affected by the small 
scale motions that are filtered out, and there is no way to recover the lost 
information; all attempts to do so failed. 

Most of the results on large eddy simulation are taken from a report by 
Mansour et al, (1978). 
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3. Anisotropic turbulence 

Anisotropic turbulence (turbulence In which the fluctuating components 

9 9 9 

are unequal so that u^ u^ U 3 ) usually returns to an Isotropic state if 

not strained In any way. However, It Is possible for the flow to become even 

2 2 

more anisotropic. [^us, If the large scales are such that u. < u» and the 

2 2 ^~~2 ^~2 
small scales have u^ > U 2 but the total field Is such that u^ < U 2 , It Is 

quite likely that the turbulence will become more anisotropic with time. This 

Is not the case In most flows, however, and the aiaumptlon that anisotropic 

turbulence tends to return to an Isotropic state Is reasonable In most flows 

of Interest. 

In the laboratory, anisotropic turbulence Is usually created by straining 
the flow and then allowing the anisotropic turbulence to relax In the absence 
of strain. The alternative approach of using the anisotropy of turbulence 
created by grids has not been successful. The apparent reason Is the one 
mentioned above — the anisotropy resides mainly In the large scales, and the 
flow may become more rather Chan less anisotropic. Creating anisotropy by 
straining an initially isotropic field distributes the anisotropy over the 
range of scales and Is thus better behaved. 

Simulations can emulate either of the above methods. One can simply 
create an Initial field In which the components of the velocity fluctuations 
are unequal, or one can strain an initially isotropic field to produce the 
anisotropy. Because one has control of the anisotropy as a function of the 
scale size in the Initial conditions In a simulation, there Is no Important 
factor favoring one method over the other. The method of creating an aniso- 
tropic initial field is preferred, as it is the simpler approach. 

Full simulations of homogeneous anisotropic turbulence were made by 
Schumann and Herring (1976) using the method suggested above. Some of their 
results are shown in Fig. 5.7. We see that their flow does Indeed relax 
toward isotropy. The tendency of the dissipation and pressure-strain terms 
toward their values in the isotropic flow is also evident. One should note, 
however, that the calculation was done on a 32^ mesh. All of the ^^uantities 
averaged in Fig. 5.7 fluctuate very strongly in both space and time, and it is 
likely that nearly all of the contribution to the mean values comes from a few 
small regions in which the fluctuations are very intense, this statement is 
based on some of the author’s unpublished work. It is therefore likely that 
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the uncertainty In the reported values is quite large. This is true in some 
of the other flows that we shall look at as well. 


Schumann and (terrlng used their results to test two versions of Rotta's 
model for the return to Isotropy. This model assumes that the pressure-strain 
term can be represented by: 


3u* 

< P’ 
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au’ 

3x, 


> = ^ • - K < u'u' > - - 6 < u'u' > 

ij ^ “i“j 1 ij k k 


(5.3) 


where one model assumes K “ Ce/q , and the other assumes K ■ C'q/L, where 
e is the dissipation rate and L is the Integral scale. The brackets < > 
represent an average over the computational field which is assumed equivalent 
to an experimental time average. As can be seen from the figure, there is 
considerable variation in the "constant" obtained from the various runs. 
Clearly, this Indicates that something may be wrong with the model. Schumann 
and Herring were not able to discern any epiisistent Reynolds number effects in 
their results. 


4. Rotating Ihrbulence 

The effects of rotation on turbulence are subtle and complex. In the 
equations of motion, the only appearance of the rotation is via the Coriolis 
force; the centrifugal force can be transformed away. One effect of the 

Coriolis force is to redistribute the kinetic energy among the components of 

the turbulence normal to the axis of rotation. The Coriolis force does not 

appear explicitly in the equation for the turbulent kinetic energy. Neverthe- 
less, rotation has a profound effect on turbulence and, especially, on its 

rate of production (cf. Ferzlger and Shaanan (1976)). In shear flows, rota- 
tion may in fact stabilize the flow, there is evidence that it can cause 

relaminarization of a turbulent boundary layer. It can also destabilize; the 

•• 

well-known Taylor-Gortler Instability is a prime example of this. 

The effect of rotation on isotropic turbulence is even more subtle. It 
seems quite likely that the principal effect is the conversion of turbulent 
energy into inertial waves— wavet that propagate principally along the axis of 
rotation and which are not dissipated except near walls. 

Experimentally, the study of the interaction of rotation and turbulence 
is very difficult. One major difficulty (about which we shall say more later) 
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is tliui the fluid must be set into rotation before it passes tbrougb a grid 
that generates tutbuience, ibis is u consequence of the Uelmboitz theorem. 
Hiree experiments have been petiormed. Ibbotson and 'I'cittun (19b7) found a 
faster decay of the turbuience wl»en the fluid was rotating, while Traugott 
(19bb) found a decrease in tl»e decoy rate. Hie latest expeirmenC, and the one 
that is generally regarded as the best, was done by Wigeland and Nagib (1978). 
Tliey found cases which went in both directions; however, the predominant 
effect was a decrease in the decay rate. 

Since the source of the effects observed in tlie experiment was unknown, 
preliminary calculations using large eddy simulation on a Ib^ grid were 
made. A aeries of simulations using the identical initial condition with 
various rotation rates was made. llie, results, sliown in Pig. 5.9, indicate 
that the predominant effect of the rotation may be to decrease the rate of 
decay of the turbulence, but there is unusual behavior, particularly at the 
early times. This is similar to the behavior observed by Wigeland and Nagib, 
but a detailed compactoon is impossible. 

On the bv'sis of tiiese results, it was surmised timt rotation decreases 
the rate of dissipation but that this effect is masked by other effects in the 
early development of the flow. In order to clieck this hypothesis, we luade 

full simulations of an experiment that is impossible to do in the laboratory. 
We allowed the turbulence to develop without rotation for a short time; this 
is identical to the initial stages of an isotropic turbulence experiment. When 
the turbulence had developed into a physically realistic field (see the pre- 
ceding section for details), the rotation was "turned on." Under these 
conditions, it was found that increasing the rotation rate always decreased 
the rate of decay of the turbulence. The results are shown in Fig. 5.10. 

It appears that the anomalous effects found in the experiments are caused 
by interactions of the rotation with the thin shear layers produced by the 
turbulence-producing grid, and similar affects can be produced in the simu- 
lation. Tltafiie are Impossible to avoid in the laboratory. in some of the 
experiments, interactions with the walls also play an important role. 

it was also possible to search for the cause of tlie effect, it was found 
that the turbulence remains nearly isotropic, so the decrease in the rate of 
dissipation must be due to an increaue in the length scales. Since the length 
scales are readily computed in these simulations, this was easily checked, and 
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It was found that there Is a large Increase in the length scale in the direc- 
tion of the rotation axis. A theoretical explanation for this (based on the 
properties of inertial waves) , was given and a modification of the model was 
offered. 

5. Strained Thrbulence 

We now come to flows in which there is turbulence production. In both 
the strained and sheared turbulence experiments, the turbulence decays for a 
short time after the start of the flow and then Increases with time. Ttie 
length scales of the turbulence Increase more rapidly than in the unstrained 
decaying isotropic flow. All of this makes these flows interesting objects of 
study. 

In the laboratory, strsltied turbulence is created by first producing iso- 
tropic turbulence with a grid, in the same way as in tlje experiments described 
earlier. Tlie turbulence is allowed to develop for a short time and is then 
made to pass through the test section. In some test sections, the cross- 
sectional area is kept constant but the aspect ratio in the plane normal to 
the flow is changed; the effect is to exert plane strain on the turbulence. In 
other experiments, the test section is a contraction, and the turbulence is 
compressed in the two directions normal to the flow and stretched in the 
streamwise direction; the result is axisymmetric strain. 

To simulate these flows numerically, an isotropic turbulent flow field is 
created in the same manner as for the previous flows. The effect of the 
strain is turned on immediately, and the flow is allowed to develop. In order 
to simulate this flow correctly, it is necessary to use a straining coordinate 
system, one which moves with the mean flow that produces the strain. This is 
necessary because one of the terms due to the applied strain does not permit 
the application of periodic boundary conditions, the transformation removes 
this term. For the details of this transformation see Rogallo (1977). 

Use of this transformation also introduces a difficulty. After some 
time, the strained coordinate system becomes quite thin in the direction which 
is being compressed. Whea this happens, the length scales in that direction 
become appreciable compared to the size of the computational domain in that 
direction. As a result, periodic boundary conditions are no longer valid, 
and the computation has to be stopped. This happens when the total strain 
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exp(St) 1 , where S ■ 3u/dx. The problem can be partially alleviated by 
starting with a coordinate system that is distorted in the other direction, 
Thua the flow contains three periods similar to those found in the flows 
described above. First there is a development period; this is followed by a 
period in which the flvw is physically realistic, finally, there is a period 
in which the simulation is invalid, and the calculation must be stopped. 

The detailed behavior of strained turbulence is dependent on the initial 
conditions. However, the trends are the same in all cases. As in the experi- 
ments, the turbulent kinetic energy decays until the turbulence becomes orga- 
nized; then the production of turbulence increases and, somewhat later, so 
does the kinetic energy of the turbulence. As can be seen in Fig. 5.11, the 
turbulence becomes highly anisotropic. The fluctuations in the direction 
being compressed (the x^-directlon for the case shown in Fig. 5.11) increase 
most rapidly, while the fluctuations in the stretched direction (X 2 ) con- 
tinue to decrease. The off-diagonal components of the Reynolds stress tensor 
are all zero in this flow. 

The results of this computation could be used to test Reynolds-avv^raged 
models, but they have not been used for this purpose. The reasons are that 
the majority of engineering flows are shear flows, and sheared homogeneous 
turbulence seems more appropriate for this purpose and that the experimental 
data can be used as well. For this reason, Reynolds-averaged models are 
deferred to the following section. 

McMillan and Ferzlger (1980) have used strained turbulence simulations 
for checking subgrid scale models. Tliey found that the Smagor insky model 
becomes less accurate as the flow is strained. The correlation between the 
exact and model results drops from the already low value of 0. 3-0.4 to nearly 
zero. However, the scale similarity model proposed in Chapter 3 is nearly 
equally valid with or without strain. 

In a few cases the correlation between the exact stress and the Smago- 
rlnsky model becomes negative. On further investigation, it is found that, if 
the strain rate is high and maintained for a long time, the energy flow is 
from the small scales to larger scales, i.e., from the unresolved or subgril 
scales to the larger or resolved scales. Uils seems to be physically correct, 
although it has not been reported in any of the experimental results of which 
we are aware. it appears that the smallest scale of the turbulence may be 
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determined by the strain rate rather ttian the viscosity* Direct evidence of a 
similar phenonemon in sheared turbulence will be presented in the following 
section. 

6. Sheared TYirbulence 

Homogeneous turbulence interacting with mean shear behaves in a manner 
very similar to strained turbulence. One can regard shear as a combination of 
strain and rotation; the effect of the rotational component is to weaken the 
effect of the strain somewhat. The behavior with time is qualitatively simi- 
lar to that for the strain case; after a period of decay, the turbulent 
kinetic energy begins to increase. The anisotropy produced is such that the 
streamwlse component of velocity has the largest fluctuations and the normal 
component has the smallest fluctuations. 

Homogeneous sheared turbulence is more difficult to create in the 
laboratory than strained turbulence. The essential reason is that, because 
shear has a rotational component, it cannot be suddenly introduced into the 
flow. It has to be created along with the turbulence. The apparatus used to 
produce this flow is an array of parallel channels whose flow resistances are 
arranged so that the velocity distribution at their exits is linear in the 
direction normal to the channel walls. In this way, a flow with a straight- 
line mean velocity profile (uniform shear) is created. With careful 
adjustment, the turbulence can be made to be approximately uniform across the 
flow. The flow is then followed down the test section, and measurements of 
the turbulence quantities are made at the midplane of the test section at a 
number of stations. 

Simulation of this flow on a computer is very similar to simulation of 
strained flow. An initial Isotropic velocity field is created In the manner 
described earlier. It Is possible to let the flow relax before the shear is 
introduced, but this is not done. For this flow it is necessary to use a 
shearing coordinate system (one that moves with the applied linear mean flow) 
in order to remove the terms that forbid the use of periodic boundary condi- 
tions. Ihe deforming coordinate system is shown in Fig. 3.12. it begins as a 
Cartesian system at t ■ 0 and deforms as shown until St ■ 1/2. Af this 
point, the computational domain is on the point of becoming too narrow in the 
normal direction to support the use of periodic boundary conditions. This 
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flow permits the "remeshlng" of the coordinate system In the manner shown in 
Fig. 5.12. The shear then causes the coordinate system to become Cartesian, 
and the cycle Is begun again. With the aid of this trick, it is possible In 
principle to carry on for as long as desired. In practice, the length scales 
in the streamwlse direction eventually become too long for the computational 
domain, and one is forced to stop on this account. Sheared turbulence thus 
passes through the same three periods as strained turbulence: development, 
realistic representation of physics, and, finally, breakdown. 

Tlie detailed behavior of the flow may depend on the Initial conditions, 
but the trends are essentially independent of how the calculation is started. 
As one can see from Fig. 5.13, the behavior of the components of the turbu- 
lence is very similar to that in the strain case. It also follows the exper- 
imental trends very well. 

McMillan and Ferziger (1980) used the results of direct simulations of 
sheared turbulence as the basis of tests of subgrid scale models. The find- 
ings differed in no important respect from thone found for strained turbu- 
lencej for this reason, we shall not give them here. However, we point out 
that the transfer of energy from the smallest scales to larger scales was 
noted in this case as well. Further evidence for this will be given below. 

Let us look at the results of the simulations in somewhat more detail. 
Many of the results are those of Feiereisen et al. (1981), Shlrani et al. 
(1981), and Rogallo (1981) which have not yet been published. Only partial 
results will be given. Three-dimensional spectra of the velocity field are 
shown in Fig. 5.14. We see that there is a very strong shift of the spectrum 
of ihe normal velocity component toward low wavenumbers or large scales. Care 
is required in dealing with the integral scales. They are the integrals of 
two-point correlation functions, some of which have regions in which they are 
negative. The negative regions can cause the integral scales to behave very 
erratically. The spectra probably show the length-scale behavior more accu- 
rately. 

The behavior of the pressure spectrum is rather remarkable. The initial 
condition has a peak at a relatively high wavenumber. The pressure spectrum 
near the end of the physically realistic period is shown in Fig. 5.15. The 
spectrum is broken into two components. The decomposition is suggested by the 
Poisson equation for the pressure; the terms on the right-hand side of that 
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equation can be claasllled according to whether or not they contain the mean 
velocity field* The component is a consequence of the applied mean 

_5 

field; it develops a k spectrum, and the peak in the spectrum moves to the 
left with advancif.ig time. The component ?2 is due to the self-interactions 
of the turbulence and is much more broad-band in nature. Tiiis has important 
consequences for pressure-strain modeling. 

Finally, wc show the time behavior of the terms that contribute to the 
spectral behavior of the turbulent kinetic energy as a function of waventtmber ; 
these are shown in Fig. 5.16. It is seen that, as expected, the production is 
mainly in the large scales or low wavenumbers, and the dissipation occurs at 
higher wavenumbers. Finally, we note that the transfer term, which redistri- 
butes energy among the wavenumbers is negative at low wavenumbers (indicating 
a transfer away from the largt:> scales) and becomes positive at higher wave- 
numbers. All of this is as anticipated. The surprise is that the transfer 
again becomes negative at the highest wavenumbers, indicating that the trans- 
fer is from both ends of the spectrum to the center. This can be taken to be 
a confirmation of the finding of McMillan and Ferzlger discussed earlier. 

Let us now look at some of the applications of these results to Keynolds- 
averaged modeling. Since this is the first application of this type in this 
report, we should first look at the possibilities. The time averages can be 
replaced by averages over the flow field. Although the number of mesh points 
is large (64 ■ 262,144), they cannot be regarded as statistically indepen- 
dent. A more realistic measure of statistical reliability is the number of 
large eddies captured in the computational domain. There are several ways to 
measure this — none of them exact — but the number of large eddies is small 
enough tliat the statistical reliability of the results is not very high. A 
good test of their validity is to compare results obtained from two simu- 
lations which are Identical except for the set of random numbers used to 
initialize them. 

From each realization of a shear turbulent flow, we may compute the 
averaged quantities as a function of time. Since the quantities vary slowly, 
the values at neighboring times are not independent, and should not be treated 
as if they are. For this reason, we chose to analyze the flow fields only at 
those times at which the grid is Cartesian. This is also convenient computa- 
tionally. The result is that we have the averaged quantities that need 
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to be modeled at three or four time steps for each of several realizations* 
The data sets are thus much smaller than those used in subgrid scale model 
testing, and the kinds of tests performed need to account for this. Further' 
more, one needs to consider the effects of changes in the basic parameters of 
the flow. 


These flows contain two independent nondimenslonal parameters. Hie first 
is the Reynolds number. There are several length scales on which a Reynolds 
number can be based. The integral scale suffers from the difficulties de- 
scribed earlier, and we have used the microscale instead. The two should be 
related (possibly as a function of Reynolds number), so it does not matter 

much which length scale Is used; however, if we try to apply the results to 

other flows, the choice of length scale may be very important. Tlie second 

nondimenslonal number is the ratio SL/q, where S is the applied mean shear 

rate and q and L are the velocity and integral length scales. We call 
this parameter the shear number, and it measures the ratio of an eddy time 
scale to the time scale Imposed on the flow. It can also be shown that the 
shear number is proportional to the ratio of production to dissipation. 


From the results of a simulation, one can compute the Reynolds shear 
stress < U|^U 2 >. This is just a single quantity, and one cannot test eddy 


Fddy viscosity models could be tested by 

I 2 

asking whether the Reynolds stress tensor, R^^*<UjU^>--jq6 


viscosity models using it alone, 
asking whether the Reynolds str( 
proportional to the rate of strain tensor 


ij 


■i“j 


ij’ 


is 


9u 9u , 

Ij “ ^'^T^ij “ 

j i 

Since the model could not be tested directly, we computed the "constant” in 
the model defined by 


< u^u^ > 
3U^/3x2 


CqL 


( 5 . 4 ) 


and correlatee it as a function of the two nondimenslonal parameters given 
above. The result showed that C is nearly inverse to the shear number, 
which is equivalent to saying that < uj^u^ q > is nearly constant, note 
that this result is incompatible with C's being a true constant. On further 
Investigation, it was found that all components of the Reynolds stress anisot- 
ropy tensor: 
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(5.5) 



< u^uj > _ i 

2 ■ 7 

q 


appear to become constant at long times in homogeneous shear flow It is 

Impossible to carry the calculation far enough to determine whether this is 
really the case or whether the b^j simply change very slowly in the later 
stages of this flow. We are of the opinion that there is asymptotic struc- 
tural similarity in this flow; this assumption has been the basis of some 
recent models. In many other shear flows ^ < Uj^U 2 >/< q > is approximately 

constant over a large part of the flow; for example, in the boundary layer 
this holds except for the region close to the wall. 


Another example of model testing with these simulations is provided by 
the pressure-strain terms. We showed earlier that the pressure can be 

considered to be composed of two parts, one arising from interaction of the 
turbulence with the imposed mean field and the other a purely turbulent 
quantity. The corresponding decomposition of the pressure-strain terras is 
made by many modelers. 

For the part of the pressure strain terms proportional to the mean strain 
(the "rapid" terms), one can show that, if one allows only terms which are 
linear in the anisotropy of the Reynolds stress, the model contains only a 
single constant, which for the 1,1 component can be written (Reynolds 

(1976)): 


< P 


3u 
1 3x 


A. 3U, 

2 1 + -5-3: < '‘i_U2 > 3 ^ 


(5.6) 


There are similar expressions for the other components. Given the computed 
values of the rapid part of the pressure-strain term, we can calculate a value 
of the "constant" for each of the four tensor components that are nonzero. If 
the model is correct, the values obtained should be the same for each tensor 
index and all realizations. The results showed that the "constant" is nearly 
independent of the Reynolds and shear numbers, but it varies by a factor of 
nearly seven among the various components of the tensor. These results show a 
deficiency in the model and suggest that an improved model should be possible, 
but we have so far been unable to suggest one. 
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The part of the pressure-atrain term that results from purely turbulent 
Interactions (the "Rotta" term) are usually modeled by: 

3u, 3u. 

^ ^2 ^ ^ ^ ‘’ij 

This model Is based on the notion that the effect of these terms Is to return 
the flow to Isotropy. It, too, Is easily tested by the method used for the 
rapid term. It was found that the "constant" displays a great deal of 
variation with Reynolds number, and many of the values were below the value of 
2 required for return of the turbulence to an Isotropic state. 

Further Investigation showed that the anisotropy of the dissipation does 
not behave as had been expected. It Is generally assumed that the dissipation 
Is Isotropic at high Reynolds numbers but may be anisotropic at low Reynolds 
numbers. Thus we expected to find a strong Reynolds number dependence of the 
anisotropy of the dissipation. In fact, we found almost no variation with 
microscale Reynolds number In the range from 10 to 100 (see Fig. 5.17). This 
does not mean that the dissipation cannot become Isotropic at still higher 
Reynolds numbers, but It does suggest that the assumption of Isotrooy may be 
questioned. 

Since the anisotropic component of the dissipation acts to reduce the 
Isotropy of the Reynolds stress tensor, it should be Included with the 
pressure-strain term. When the combined terms are modeled. It is found that 
the variation of the "constant" with Reynolds number is greatly reduced (see 
Fig. 5.18), and the model Is fairly good. Modelers who assumed the dissipa- 
tion to be isotropic have gotten reasonably good results because the aniso- 
tropy of the dissipation is implicitly Included in their models. 

This Is a sample of some of the results obtained by Felerelsen et al. 
(1981) and Shiranl et al. (1981). The reader is referred to those reports and 
forthcoming papers for more complete details. 

7 . Compressible 'Ibrbulence 

It Is possible to make a compressible version of the homogeneous turbu- 
lent shear flow treated in the preceding section. One need only make the 
velocity gradient large enough that the velocity difference across a large 
eddy is a significant fraction of the sound speed. It is not possible to 
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produce this flow in the laboratory, the large velocity differences would make 
it impossible to maintain homogeneity. Tltis is unfortunate, b< 3 cause it means 
that we have to believe the results of the calculation Wi.tl’out experimental 
verification. We can, however, check the results at low Mach number against 
the incompressible experiments. 

To compute this flow, the major change We need to make from the incom- 
pressible case is that the full set of compressible equations must be used. 
One can show that a linear velocity profile is a solution to the steady equa- 
tions, and this solution can serve as the source of the shear imposed on the 
turbulence. In compressible computations (cf. Uallhaus (19B0)), it is custom- 
ary to use the continuity, momentum, and energy equations in conservation 
form; the dependent variable in the energy equations is usually the total 
energy (stagnation enthalpy). However, in the present case, this equation 
cannot be used without destroying the homogeneity (Feiereisen et al. (19B1)), 
80 we are forced to treat the enthalpy as one of the primary dependent vari- 
ables. 

The most popular numerical methods for the compressible equations are 
designed to relax the solution to a steady state as quickly as possible. Tliey 
are not time-accurate; that Is, they do not produce an accurate picture of the 
relaxation to steady state, and therefore they cannot be used for the purpose 
we have in mind. Instead, we have used a standard explicit method. The 
fourth-order Runge-Kutta method was chosen. The fact that all of the compres- 
sible equations contain time derivatives means that one does not need to solve 
a special equation for the pressure. All variables are advanced in time, the 
variables which are not explicitly computed from the differential equations 
are obtained from equations of state. 

Morkovin (1963) hypothesized that compressible turbulence behaves very 
much like incompressible turbulence, and most models are based on this assump- 
tion. For most of the quantities in homogeneous turbulent shear flow, this 
hypothesis turns out to be correct. Most of the differences between the two 
cases are small, so we shall concentrate on the few cases in which the differ- 
ences are significant. 

The major difference between the incompressible and compressible flows 
(at least when the turbulence Mach number is not too large) is due to the 
appearance of acoustic waves in the latter case. The acoustic waves that are 
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moBt apparent are thoae propagating normal to the shear, and we expect the 
quantities which can be affected by acoustic waves to show the most important 
differences from the incompressible case. The largest ctmnge is in the 
fluctuating velocity component normal to the shear; it is reduced relative to 
the incompressible case. 

The most striking difference between the two flows is in the pressure and 
the terms associated with it. In the incompressible case, the pressure was 
decomposed into two parts: one arising from the mean flow that produces the 
shear and another tliat is entirely due to the turbulence. In the compressible 
case, there is a third term due to the presence of acoustic waves. This term 
turns out to be significant even at fairly low Mach numbers. 

Of course, the pressure-strain terms are also affected in the same way; 
there are now three of them. It turns out that che third term behaves like 
the rapid term — 'the one due to the mean shear — and can therefore be combined 
with it. However, the "constant” is now a function of the turbulent Mach 
number in addition to the two dimensionless parameters of the incompressible 
case — the Reynolds and shear numbers. The resulting constant was fit as a 
function of these three parameters. The results are shown in Fig. 5.19, and 
the Mach number dependence is found to be significant. 

Further details and results for this flow can be found in the report of 
Feiereisen et al. (1981). 

8. Mixing of a Passive Scalar 

By definition, a passive scalar is any quantity that can be convected by 
a flow and diffuse through it without affecting the velocity field. There are 
many applications that require knowledge of how a passive scalar behaves, any 
problem in which heat or mass transfer is important is in this class. Under- 
standing the mixing of a passive scalar is also a preliminary to handling 
reacting flows, including combustion. 

A passive scalar could be introduced into any flow treated in this chap- 
ter. In fact, only two of these have been done experimentally, these are 
isotropic turbulence and sheared homogeneous turbulence, so these are the 
cases which have been simulated. One also has to decide whether the scalar 
has a mean component or not. In the experiments, isotropic turbulence has 
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been meaeured without: a naan gradient of tlie acaiar , and the shear flow haa 
been performed both with and without a mean «?alar gradient « To facilitate 
compariaon with theae experimenta, iaotropic turbulence was simulated with an 
isotropic scalar field, and the shear flow had a mean scalar gradient. 

The equation describing the scalar concentration la : 


3C 

3t 


+ 




a^c 


(5.8) 


If there is a mean scalar field, it is subtracted from the total scalar field 
to obtain an equation for the fluctuating scalar field. The velocity field is 
also decomposed into its mean and fluctuating parts. The resulting equation 
for the scalar fluctuations has the same difficulty as the equation for the 
velocity fleld-'-'the mean shear and mean scalar gradient terms do not admit the 
use of periodic boundary conditions. To remove this problem, the coordinate 
traneformatlor made for the momentum equations has to be made here as well. 
It is possible to compute the velocity field prior to the computation of the 
scalar field, but this would require storing an enormous data set on tape and 
transferring it back into the machine as needed. For this reason, the 
velocity and scalar fields were computed simultaneously. The numerical 
methods used for the scalar field are identical to those used for the velocity 
field. 


In the case of the isotropic field, the most important items to study are 
the decay rates of the velocity and scalar fields. The scalar field follows a 
decay law similar to Eq . (5.2): 


c'^ - B(t-t^)"“ 


(5.9) 


where c' is the fluctuating part of the scalar field, i.e., c*<C>+c'. 
We wish to look at the ratio m/n. The parameters on which this ratio depends 
are the Reynolds number and the Prandtl or Schmidt number, which is the ratio 
of kinematic viscosity to dlffusivlty (Sc “ v/D) . It was found that the 
scalar decays more rapidly than the velocity field when the Schmidt number is 
less chan unity and more slowly than the velocity field when the Schmidt num- 
ber is greater than unity^ this is no surprise. The dependence of the ratio 
m/n on Reynolds number also depends on whether the Schmidt number is less 
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than or greater than unity. For Sc < 1 , it is found that the ratio m/n 
decreases with increasing Re, and vice versa for Sc > 1 . 

The cases which include shear and a mean gradient of the scalar were ana- 
lyzed in a manner similar to that used for the homogeneous shear flow. An 
Important result is that the behavior of the scalar field becomes Independent 
of the initial conditions after a short time. Its properties depend almost 
entirely on the velocity field and the mean gradient of the scalar. 

The next quantity studied was the scalar fltut < u^^c >. This quantity is 
usually modeled by gradient diffusion: 


< u^c > 


a < C > 


ij 


J 


( 5 . 10 ) 


In the standard case,, the concentration gradient is in the same direction as 
the velocity gradient; the nonzero gradients are 3Uj^/3x2 3C/3x2 and 

there are two nonzero eddy dlf f uslvlties , Dj_ 2 and D 22 ’ The important one 
in most applications is ^22’ computed for a number of different val- 

ues of the dimensionless parameters of the flow. One can form the turbulent 
Prandtl/Schmldt number, by taking the ratio of the eddy viscosity to 

the eddy dlffusivity. A number of models have been proposed for Pr^, and 

the ones that were recommended most strongly in the literature were tested. 
None of them was found to be very accurate. A new model was constructed which 
gives in terms of j » the anisotropy of the Reynolds stress tensor. 

Although this model models a low-order quantity in terms of a higher-order 
quantity, it can be made into a useful correlation by using other correla- 
tionsi this model was found to be a significant improvement over the ones 

suggested in the literature. A test of the new correlation is shown in Fig. 

5 . 20 . 


It is also possible to compute the other nonzero elements of 
to the design of the computer program, this was not done for the full range of 
cases for which D 22 was computed. Also, since the elements of the diffusiv- 
Ity tensor depend on the nondimenslonal parameters and, because it was not 
possible to match the Reynolds number used in the experiments, a quantitative 
comparison with experiment is not possible,. However, the results are in good 
qualitative agreement with the data. 
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We also correlated the mean-square scalar fluctuations as a function of 
the nondimenslonal parameters. Ttie principal finding was that, to a good 
approximation , 


< c > 

■5OTT 




(5.11) 


One can also construct models for the scalar field based on the ideas 
used for the velocity field. In particular, one can derive equations for 
< c^ > and < CU2 > , which are similar to the Reynolds stress equations. 
The terms in them that are most difficult to model are the correlations be- 
tween the fluctuating pressure field and the gradient of the fluctuating 
concentration < p 3c/ 3x^ > . They are analogous to the pressure-strain 
terms, and models for them can be based on models used for the latter. In 
particular, the pressure decomposition used in deriving pressure-strain models 
can be used here as well . 


The model for the rapid term (the one containing the pressure derived 
from the mean shear) contains no adjustable constants. However, we introduced 
an arbitrary multiplicative constant and found good agreement between the 
exact and model results. The constant was found to be approximately 0.5, 
Indicating that the arguments made in deriving the model are deficient. 
Another model suggested by Lumley to overcome some of the undesirable prop- 
erties of the first model was tested and found to be less accurate than the 
first model. 


The term arising from the component of the pressure that depends entirely 
on the turbulence was modeled by an analog to Eq . (5.7). The results show 
this model to be quite good — ^about as good as the modified Rotta model de- 
scribed earlier. 
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Chapter VI 
FREE SHEAR FLlr*! 


I . Overview 

Free shear Clows are one of the classes of flows of major technological 
interest. They occur in many kinds of devices, and we shall begin this 
chapter by briefly describing the types of free shear flows. 

Free shear flows can be divided into three major categories, there are 
also more complex cases. The thr»^e major types are: 

^ • Mixing layer . This is tlie flow that occurs when two parallel flows 
of different velocity are brought together. In the laboratory this flow is 
created by having tite fluid of different speeds on opposite sides of a divid- 
ing plate. At the end of the plate, the two streams come into contact, and 
the thickness of the layer in which the velocity gradient occurs grows with 
downstream distance. 

2. Jet . A stream of high-velocity fluid issuing from an opening is 
called a jet. As the high-speed fluid mixes with the surrounding lower-speed 
fluid, the maximum velocity of the jet decreases, and the rate of growth of 
its thickness also decreases. The most commonly studied jets are the plane 
and round ones, but others, such as the rectangular jet, have been studied. 

3. Wake- , A wake is similsr to a jet, but it is a velocity defect in an 
otherwise uniform stream. Like the jet, the wake has decreasing velocity 
gradients with downstream distance. Most wakes result from flows around 
bodies. The wakes form by merging of the boundary layers behind the body or 
from separation of the boundary layers. 

We should also mention : 

Complex shear layers . This is not a single type of flow, but a category 
containing the flows that do not fall into the above categories. Curved jets 
and wakes are quite common. Another important flow is one in which a laminar 
boundary layer separates, creating a free shear layer. The free shear layer 
then undergoes transition to a turbulent free shear layer which grows so 
rapidly that it soon reattaches to the surface. This is a common mechanism of 
transition from a laminar to a turbulent boundry layer. 
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It is also important to dlstingulah the early pliases of free shear flows 
from the far-downstream flows. The early stages are sensitive to the initial 
conditions. Fully developed free shear layers are usually self-similar in 
nearly all of the measured variables and grow according to a power of the 
downstream distance. A majority of free shear layers occurring in applica- 
tions are of the early type ^ but fully developed cases are also of importance. 

To date, there have been large eddy and full simulations of mixing layers 
and full simulations of wakes. The jet has not yet been simulated (although 
it probably will be in the near future) . Complex free shear flows liave also 
not yet been attempted. Thus we shall devote the rest of this ctiapter to the 
mixing layer and the wake. 

Nearly all laboratory free shear flows develop with downstream distance. 
It is much easier to simulate a layer that develops in time. One must be very 
careful in comparing the two cases . Consider the mixing layer . Fluid ele- 
ments on the two sides of the laboratory shear layer have been in the flow for 
differing amounts of time. As a result, the development of the flow is not 
symmetric and the plane on which the mean velocity is the average of the two 
free stream velocities is inclined. The simulated mixing layer is, however, 
symmetric. The two flows may be compared if, in the laboratory flow, the 
velocity difference across the flow is small compared to the average veloc- 
ity. This experiment requires a long apparatus, but cases exist which meet 
this criterion fairly well, and these are the ones to which the simulations 
should be compared. 

2 . Mixing Layer 

As discussed above, this is the simplest of all of the free shear flows. 
Despite the apparent simplicity of this flow and the large number of experi- 
ments that have measured it, there is still controversy about it. Let us 
consider the fully developed mixing layer first and the t\ ?’asitional case 
later . 

It is generally agreed that the velocity profile of the fully developed 
mixing layer is self-similar, and so are the components of the Reynolds stress 
tensor. Another point of general agreement is that the growth of the free 
shear layer is linear in the fully developed region. The major point of 
disagreement in this regime of the flow concei ns the rate of growth of the 


69 


layer. For the mixing layer sketched In Fig. 6.1, the growth rate parameter 
is conventionally defined as : 


d6 „ ~ “l 

dfx " + “2 


( 6 . 1 ) 


The measured values of 0 cover the range 0.0U-U.16, a much wider range 
than would be expected for a flow this simple. Birch (1980) recently reviewed 
the data and believes that there Is a single correct value of this parameter, 
which he believes to be 0.113. iiowever , no reason was given for the spread in 
the data. 

There is less agreement about the early stages of the shear layer. One 
group, including Roshko and his coworkers and Browand and his coworkers, among 
others, believes that this part of the flow is essentially two-dimensional. 
In this view, the initial laminar shear layer rolls up into two-dimensional 
vortices, which then agglomerate or pair to form larger vortices of the same 
type (with larger spacing) . This process has been observed to continue for 
several pairings. At this point the flow reaches Che end of the apparatus. 
According to this view, the important process in the growth of the mixing 
layer is the pairing of the vortices., However, there is evidence that stream- 
wise vortices form in this flow. This is a kind of three-dimensionality, but 
it is quite regular rather than chaotic. 

The other view, held by Bradshaw and others , is that the mixing layer is 
normally strongly three-dlmenslonul and chaotic. According to this picture, 
the highly two-dimensional layers that some experimenters have observed are 
the result of careful arrangement of the initial conditions and design of the 
experimental apparatus. 

Large eddy simulations of the mixing layer were made by Mansour et al . 
(1978). They used the vortlcity equations rather than the primitive equations 
because the vortlcity is confined to a relatively narrow region of the flow. 
In fact, it appears that it makes little difference which set of equations 
is used. The subgrid scale model had Co be modified to account for this 
change. At the top and bottom of the computational region, no stress boundary 
conditions (see Section 4.5) were applied. Fourier sine and cosine transforms 
were used in the normal direction. 
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This work showed that it is possible to explain the rapid growth of the 
mixing layer by vortex pairing. The flow was begun with an initial condition 
that contained well defined two-dimensional vortices. Although there were 
only two vortices in the computational domain^ the boundary conditions imply 
that they are part of an Infinite array. Various perturbations to a regular 
vortex array were tried. It was found that small perturbations would cause 
the vortices to pair. Naturally, the pairing occurred more rapidly when the 
perturbation was larger. Surprisingly, it was found that the mean veloclt)*' 
profile (defined by averaging the velocity over a plane) was self-similar and 
that the growth of the momentum thickness of the layer was very nearly linear. 

A number of three-dimensional perturbations on this basic flow were also 
made. First, small, random, three-dimensional disturbances were added to the 
Initial conditions. The three-dimensionality was somewhat amplified by the 
pairing process, but there were only minor changes in the overall properties 
of the flow. Another variation was produced by the addition of streamwlse 
vortices to the initial condition. The streamwlse vortices were distorted in 
the pairing process , and they produced slight kinks in the large two- 
dimensional vortices that result from the pairing. It was conjectured that 
the kinks would produce larger-scale instability of the mixing layer and would 
then lead to considerable three-dimensionality, but this could not be demon- 
strated because the number of grid points was severely limited, 

A simulation of the initial stages of the mixing layer was made by Cain 
et al. (1981). This simulation used numerical methods described in Section 
4.5. The transformation of an infinite region to a finite one was used, and 
the modified Fourier method of taking spatial derivatives in the normal direc- 
tion was used. The initial profile was a laminar mixing layer with a small 
random disturbance; the disturbance was strongest on the center plane of the 
layer . 

The results turned out well. Use of the coordinate transformation and 
the Fourier method allowed the method to be continued until the original layer 
had grown by nearly a factor of ten in some cases. No effect of image layers 
was found, and, in most cases, the calculation was stopped only because the 
layer developed horizontal scales which were too large to permit application 
of periodic boundary conditions. 
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Several variations In the computational method were tried. A full 
simulation was made; this calculation was stopped because the turbulence in 
the small scales became too strong. The calculation was repeated with fil- 
tering, but no subgrid scale model; the problem with the small scales disap- 
peared, and the calculations could be carried almost twice as far in time, at 
which point the difficulty with the large scales appeared. A final calcula- 
tion with both filtering and the subgrid scale model was made; it differed 
only a little from the preceding case. Thus, most of the results were ob- 
tained with filtering but no model. 

Simulations were made with three levels of initial disturbance. In the 
low initial turbulence cases, the turbulence intensity was four orders of 
magnitude smaller than that of a fully developed turbulent layer, this might 
represent the behavior of a mixing layer produced from laminar boundary 
layers. The medium initial turbulence level was two orders of magnitude 
stronger. The high initial turbulence level cases started with turbulence 
intensities nearly those of the fully developed layer; these might represent a 
mixing layer produced from turbulent boundary layers. Cases which differed 
only in the set of random numbers used to generate the initial conditions were 
also run. 

The results show chat the low-turbulence cases produced a layer in which 
the momentum thickness grew very slowly at first but, after a latency period, 
grew linearly with time at a rate similar to that observed in experiments. 
The medium-level case gave a shorter latency period and a slightly slower rate 
of growth at later times. Finally, the high- turbulence level cases gave 
almost no latency period at all but a still slower growth rate. These results 
are in qualitative agreement with experimental data. They are shown in Fig. 
6.2, 6.3, and 6 .4 . 

All of the cases have mean velocity profiles that are seif-similar. 

The growth of the turbulence level on the center plane of the mixing 
layer is shown in Figs. 6.5, 6.6, and 6.7. In the low-turbulence case, the 
turbulent kinetic energy grows exponentially in the early stages and then 
levels off; the exponential growth rate is close to that of the most rapidly 
growing mode according to linear stability theory. The medium initial turbu- 
lence cases show similar growth, but the exponential period does not last as 
long. The high initial turbulence cases grow only slowly as they begin near 
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the level for a fully developed layer. In all of the cases the kinetic energy 
of tJie turbulence overshoots the value for the fully developed layer before 
settling down. This has been observed in some experiments. 

The profile of the turbulent kinetic energy is shown for a typical case 
in Fig. 6.8. The initial profile is too broad compared to the fully developed 
profile. This is corrected^ but the profile becomes too Chin before the final 
state is reached. 

The simulations were also used as Che basis for flow visualix.ations . A 
grid of “dye lines" was placed on Cte center plane of the flow at the initial 
time. The ones in the streamwise direction are essentially vortex lines in 
the low-intensity cases and remain so by Helmholtz's theorem. The dye lines 
are moved with the flow, and pictures ace drawn at various times. The initial 
picture is shown in Fig. 6.9, and the final result is shown for two different 
initial fields in Figs. 6.10 and 6.11. It is clear that the layer has rolled 
up into vortices, but they are much more two-dimensional in one case than the 
other. We believe Chat the three-dimensional shear layer does roll up into 
vortical structures, but that these structures do not have spanwise uniformity 
except when precautions are taken to insure chat the three-dimensional distur- 
bances are weaker than the two-dimensional ones. 

The above results were taken from the report of Cain ec al. (1981). 

Two-dimensional simulations of the mixing layer were made by Patnaik, 
Sherman, and Corcos (1976), Acton (1976), Knight (1375), Ashurst (1979), and 
Riley and Metcalfe (1980), among others. In these simulations, the shea.r 
layer rolls up into an array of vortices. The principal object of these 
studies was the determination of the effect of initial perturbations on the 
speed and nature of the rollup of Che layer. These papers contain interesting 
results, but, as they ace essentially outside the topic of this report, they 
are not covered here in detail. 

Full simulations of the turbulent mixing layer were made by Riley and 
Metcalfe (1980a ,b). These simulations are similar to the work of Cain et al . , 
which they predated. Their calculations are performed at low Reynolds number 
so that no subgrid scale model is required. Their initial condition is eirai- 
lar to the high initial energy condition of Cain et al., but they also ran a 
number of cases in which a deterministic perturbation was added to the Initial 


73 


conditions; this perturbation was the most unstable wave of linear theory. 
They observed that the layer tended to roll up Into vortices and found linear 
growth of the thickness of the layer ^ self-slmllarlty of the velocity profile, 
and, In the case with the largest number of mesh points, constancy of the 
turbulent energy In the center plane of the layer. All of these observations 
are In agreement with experiment and the computations described above. An 
Important contribution of this work Is the demonstration that the properties 
of the mixing layer can be reproduced In a simulation which contains no large 
vortical structures In the initial conditions. They also showed that the 
addition of the perturbation corresponding to the m-.ist unstable wave of linear 
theory to the Initial condition reduced the rate of growth of the layer . 

3 . Wakes 

As stated in the introduction to this chapter, wakes are flows In which 
there is a defect In the velocity profile. As a i/ake develops, the velocity 
profile widens and the velocity gradients decrease. These factors and the 
fact chat Che rate of growth of Che length scales is not as rapid in wakes as 
in mixing layers make wakes a little easier Co simulate than mixing layers. 

There are several types of wakes. The classification plays some role in 
determining how the flow will be simulated. A selt-propelled body (one that 
drives itself through the fluid) leaves a wake in which the net momentum is 
zero; the momentum added by the propulsion just equals Chat due to the drag of 
the body. On the other hand, the wake of a towed bod.y (or a body in a wind 
tunnel) has a net momentum deficit. Finally, both types of wakes can occur in 
plane, axisymmetric , and other geometric arrangements. 

The first full simulation of a momentumless wake was made by Orszag and 
Pao (1974). Their work has been extended to the simulation of towed wakes by 
Riley and Metcalfe, in a series of papers. They concentrated mainly on the 
axisymmetric wake, because most of the experimental data is for this case. 
Despite the axisymmetry of Che flow, they used a rectangular grid in their 
simulations; the axisymmetry is inserted via the initial condlcions. 

In some respects, their simulations behave very much like the simulations 
of the preceding section. As in all other flows, a short time period is 
required for the initial condition to develop into a physically realistic 
flow. During this period there is relatively little broadening of the wake 
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aad some decay of the turbulence. The higher-order atatiatlca also change 
from their original valuea during this period; in particular, the velocity- 
derivative skewness increases. Finally, tlie vorticity tendu; to concentrate 
during this phase . 

Figure 6.12 shows the decay of the maximum mean velocity and the maximum 
axial component of the turbulence. Several experiments have shown that these 
quantities decay as with downstream distance. Since the simulated 
wakes are temporally developing, the analogous behavior would have these quan- 
tities decay as The figure shows that the maximum mean velocity fol- 
lows the expected similarity behavior quite well. The turbulence decays a 
little more slowly than expected. Two different reall;£aCions of this flow are 
s hown . 

Similarity arguments suggest that the radii of the wake and of the turbu- 

1 /3 

lence profile should increase as t ' , the spatially decaying wake radius 
Increases in radius as x^^^. Figure 6.13 shows that the simulation repro- 
duces this behavior quite well. The decay of the integrated mean and turbu- 
lent energies are also well predicted. 

The velocity profiles behave in a self-similar manner after the initial 
period. They agree quite well with the measured profiles except In the wings 
of the profile; the results are shown In Fig. 6.14. The Reynolds shear stress 
is also reasonably well predicted, as are some of the higher-order statistics. 

To conclude this chapter, we note that full simulations seem to be able 
to predict free shear flows quite well. The major stumbling block to continu- 
ing the simulations further in time is the growth of the length scales with 
downstram distance or time. This can be partially cured by doing Che simula- 
tions with larger numbers of grid points. It would be more efficient Co re- 
scale the problem after some time, but no way has yet been found to do this 
without invoking very serious approximations. 
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Chapter VII 


WALL-BOUNDED FLOWS 


1 . Overview 

The last group of flows that we shall consider in detail in this report 
is the wall-bounded flows. This is the most studied single class of flows 
because of its many Important technological applications. Despite the enor- 
mous amount of analytical and experimental attention lavished on these flows, 
there remains a great deal to be done. 

The most important single flow in technological applications is the tur- 
bulent boundary layer. The standard case for this flow is the boundary layer 
in the absence of "extra rates of 8traln"--no pressure gradient, curvatur*' » 
rotation, suction, blowing, or roughness, etc. A great deal is known about 
this flow. In particular, the mean velocity profile has been well measured, 
and one can "predict" its behavior. (Quotes are used because all of the pres- 
ent prediction schemes rely heavily on experimental data and should be called 
"postdictlve" methods.) However, the mechanism by which momentum is trans- 
ferred to the wall is only partially understood. Furthermore, the Information 
that is available about the mechanism has not been used in model construction. 
Thus there is still much to do. It is hopefl that higher-level simulations can 
play a role in this, but it is clear at the outset that the task is not easy. 

It is known that the inechanisiri of momentum transfer to the wall in the 
boundary layer is connected with the flow structure observed close to the 
wall. In the near-wall region, the flow consists of alternating "streaks" of 
high- and low-speed fluid; the streaks are very long in the streamwise dlrec- 
tl ■;< and thin in the spanwise direction. Their dimensions are believed to 
scale with the shear stress, which is nearly constant in the vicinity of the 
wall; however, their size relative to the boundary layer thickness is quite 
Reynolds number-dependent. The mechanism of momentum transfer involves lift- 
ing of the low-speed streaks from the wall. When they are lifted, they are 
carried a considerable distance into the boundary layer and exchange momentum 
with the fluid they encounter there. The existence of streaks and their 
importance in the flow plays a very Important cole in the simulation of wall- 
bounded flows . 
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The boundary layer is made up of at least three sublayers. Tliete is an 
inner layer in which the viscosity plays an important role (the viscous sub" 
layer), here the length scales are dependent on the shear stress and are small 
compared to the boundary layer thickness. The outer region of the flow is 
essentially inviscid and behaves much like a free shear low. In fact, it is 
frequently called the *'wake" region, in the wake region the length scales are 
approximately 0,1 of the boundary layer thickness. Between these two re- 
gions is one in which the shear stress is nearly constant and the viscosity is 
not important. In this region, the mean velocity has a logarithmic profile, 
and it is called the logaritlanic or buffer region, here the length scales 
increase linearly with distance from the wall. Tliis knowledge is very impor- 
tant in higher-level simulations of these flows. 

The turbulent boundary layer increases in size with downstream distance. 
This is difficult for higher-level simulations to handle at the present time. 
One can consider a temporally developing boundary layer, this has been done 
and will be described in the last section of this chapter. Unfortunately, the 
velocity profile of the time-developing layer is different from that of the 
spatially developing layer, and the difference is significant because wall- 
bounded flows are quite sensitive to small changes in the velocity profile. 

Most of the attention to date has been given to turbulent channel flow. 
It is the ideal choice for simulation, because it is the one true "equilib- 
rium" flow of the class. It reaches a state at which none of its properties 
changes with do'wnstream distance. Despite this, the physics of the near-wail 
flow Is similar to that of Che boundary layer. 'Dius this flow can be simula- 
ted with periodic boundary conditions without making any important assumptions 
that might affect the results. Of course, one must be careful that the usvjal 
criteria needed for the application of periodic boundary conditions be main- 
tained. This flow has been simulated a number of times and will occupy the 
major part of this chapter. 

Another very important issue in wall-bounded Clows is that of transi- 
tion. Laminar boundary layers are much more stable than are laminar free 
shear flows, but transition takes place when the Reynolds number is high 
enough. Tlie ability Co delay transition would enable us to reduce the drag on 
bodies, with obvious and important consequences. This is one of the major 
reasons why transition has received so much attention. 
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Transition In boundary layera la sensitive to relatively small changes In 
the velocity profile. 'Keeping the disturbance level small can delay 
transition for a long way. On the other hand, minor disturbances, such as a 
bit of roughness, can trigger transition. 

Theory predicts that laminar channel flow is stable with respect to small 
disturbances at Reynolds numbers below about 5700. One can also show that It 
Is more unstable with respect to large disturbances, but the predicted Rey- 
nolds number of transition is smaller than the Reynolds number at which 
transition is observed to occur. An explanation this phenomenon will be 
given later in this chapter. 

The next section will cake up the computation of fully developed channel 
flow. There are two approaches to doing this, and we shall discuss them and 
give results obtained by both approaches. In particular, we shall describe 
recent results that promise to provide a great deal of interesting information 
about this flow. 

The last section of this cnapter will consider transition in wall-bounded 
flows. This problem has been done recently for both the channel and the time- 
developing boundary layer* A number of interesting results have been pro- 
duced, and there is considerable hope that still more will be forthcoming in 
the near future. 

2. Fully Developed Channel Flow 

The dynamical behavior of fully developed channel flow is similar in many 
respects to that of the boundary layer. In particular, the inner layers of 
the two flows are quite similar. 'Die major differences are that the channel 
flow requires a pressure gradient to overcome the frictional forces and that 
the channel flow has no region in which the flow is not completely turbulent, 
outside the boundary layer, the flow is potential. 

Of particular importance in the simulation of the channel flow is the 
behavior of the length scales. What makes these flows especially hard to 
simulate is the fact that the spanwise length scales are much smaller near the 
wall than in the central portion of the flow. Tlii s means that a grid that is 
well adapted to capturing the streaks near the wall will be much finer than 
necessary near the center. On the other hand, a grid which is scaled for the 
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central region vlll not be able to see the ntreake at all. Ihe viirlation in 
the length scales in the direction normal to the wall is less serious , because 
a variable grid size can be used in this direction. 

TVo approaches have been taken to simulate channel tlow. in the first 
method^ which was developed by Deardorff and extended by Schumann and co- 
workers, the wall is not treated explicitly. TSile avoids much ot the dlfii- 
culty with the small-scale structures that occur near the wall, and reduces 
the amount of computation considerably* The limit of the computational domain 
is placed in the logarithmic region of the flow, because this is probably the 
best understood part of the flow. Another argument put forward for this 
method is that viscous effects prohibit the existence of an inertial subrange 
in the inner layers, but one exists in the buffer and wake regions. Tlie 
difficulty with this method is that the boundary conditions at the top and 
bottom of the computational domain are not well defined, and assumptions must 
be made. Also, this approach does not simulate much of the physics of the 
flow and cannot be used to study its structure and modeling. 


Deardorff assumed that the derivative of the streamwise velocity in the 
normal dicectin could be written as a sum of two parts, the first guarantees 
the existence of a logarithmic region, and the second is responsible for the 
fluctuations. His expression is; 


a2“ 
a u, 
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(7.1) 


where K is the von Karman constant (0.41) and Ax 2 is the distance of the 
first mesh point from the wall. This boundary condition assumes that the 
fluctuations of the velocity are the same in the normal and spanwise direc- 
tions. The validity of this assumption is open to question, Che reason that 
Deardorff gave for favoring it is that it produced reasonable results. Kim 
(private communication) has tested this boundary condition and finds it is not 
good at all. 


Schumann's assumption is that the shear stress and the velocity are In 
phase at the first mesh point; according to Kim, this assumption is also 
Inaccurate. Mathematically, his assumption is: 


< T > 
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where < u > Is the mean velocity at the first mesh point, < > is the 
mean wall shear stress, and u is the instantaneous velocity. 

For the subgrid scale model, Deardorff used the Smagorinsky model. The 
only modification that he found necessary was the reduction of the magnitude 
of the constant in the model from uhe value obtained from theory ''r the iso- 
tropic decay simulations. 

Schumann modified the model. He assumed that the subgrid scale model 
should be composes of two parts. The first is proportional to the time-mean 
velocity gradien- at the particular distance from the wall, the second is 
proportional to the deviation of the instantaneous velocity from the time- 
mean. He called these the inhomogeneous and locally isotropic components of 
the subgrid scale stress. He also used an equation for the subgrid scale 
turbulent kinetic lergy , but found that it gave no significant improvement 
over an algebraic eddy viscosity model. 

For the mean velocity profile, Schumann obtained very good results. The 
results for the components of the Reynolds stress are also quite good. Schu- 
mann also used his results for testing the Rotta model for the pressure-strain 
term. These results are shown in Figs. 7<1 and 7.2. It is interesting to 
note that the "constant" is different for the various components. However, 
one should be cautious about accepting these rsults, because the pressure is 
very sensitive to changes in the way in which the flow is computed, and we 
believe that large uncertainties must be assigned to these results. In fact, 
the results near the boundary seem to be due to the boundary conditions used. 
We shall have more to say about this below. 

Moin et al. (1978) made the first attempt to solve the channel flow prob- 
lem while treating the wall boundary conditions exactly. Doing this means 
that a nonuniform grid has to be used in the direction normal to the wall, the 
use of Chebychev polynomials is an alternative. 

One of the major difficulties with this method is that the length scales 
of the flow become very small near the wall, the local turbulence Reynolds 
number also becomes very small, and it is not clear that the Smagorinsky model 
can be used any longer. In fact, it is possible that the overall length 
scales of the turbulence will be smaller than the size of the grid in this 
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region. It Is then Improper to use the grid or filter size in the subgrid 
scale model. Instead, Muin et al. used the minimum of the L^rundtl mixing 
length and the grid size. 'Rjis modification is arbitrary but is a simple 

method that appears to work. 

Another difficulty is that the smallness of the grid tends to make 
numerical methods unstable. lliere are two nondimensionai numbers that deter- 
tilMe the stability of a numerical method* They are the Courant number 
U 2 Ax 2 /At and the viscous parameter VAt/Ax 2 * Roughly speaking, stability 
requires that both of these numbers be smaller than some constant of the order 
of unity. It turns ^ut that the viscous condition is more stringent near the 
wall, and if an explicit method were used, it would be necessary to use an 
extremely small time step. Consequently, a method which treats the viscous 

terras containing derivatives with respect to the normal coordinate implicitly 
was devised and used. Doing this meant that the normal method of solving for 
the pressure via the Poisson equation had to be abandoned. We shall briefly 
describe the revised numerical method. 

Most of the terras in the momentum equations are time-differenced using 
tlie second-order Adams-Bashf or th explicit method. The exceptions are the 
pressure gradients and tlie viscous terms containing derivatives with respect 
to the normal coordinate, which are treated by the implicit Crank-Wicolson 
method. 1V»e continuity equation, which contains no time derivatives, is 
evaluated at the nr w time step. The resulting set of equations if Fourier- 
transformed in both horizontal directions to produce a set of equations which 
are essentially ordinary differential equations with respect to the normal 
coordinate. These are finite-dif ferenced by a second-order method, and Che 
resulting set of equations is block-trldiagonal with 4x4 blocks. This 
system Is easily solved by a standard block-tr idiagonal algorithm, and, when 
the resulting functions are inverse Fourier-transformed, we have the velocity 
and pressure fields at the new time step. Kira and Moin (1979) made improve- 
ments on this method. 

The initial conditions were described in Chapter 3; t',;.y consist of a 
mean profile, solutions obtained from stability theory, and random fluctua- 
tions, TTie program was required to run for some time to ascertain Chat tlie 
turbulence wculd not decay and to develop the proper statistics. When this 
was done, it was found that the resulting velocity field contained many of the 
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features observed in the laboratory. In particular, the mean velocity profile 
was very close to the experimental one, and the fluctuating components were 
also quite close to the experimental ones. The most interesting observation 
about the results was the tendency of the fluid near the boundaries to form 
high- and low-speed streaks and for the Reynolds stress to be highly inter- 
mittent in both time and space. All of this suggests that much of the physics 
Is captured. However, the grid was not fine enough to resolve the small 

structures near the wall adequately (the "streaks" are too wide), and the 
quantitative results have to be treated cautiously. Moin et al. (1978) showed 
that this approach to simulating wall-bounded flows can succeed and indicates 
that better resolution would probably produce still better results. The 
pressure-strain correlations calculated by Moin et al. (1978) differ consid- 
erably from those of Schumann (1973). One should be very careful about 
accepting any of these results without further confirmation. Tlie preseure- 
straiu results are very sensitive to small changes in the flow, we believe 
that the trends (and the "splat" effect in particular) are correct, but the 
quantitative values are somewhat uncertain. 

Over the last three years, Kim and Moin have improved the channel flow 
calculation in a number of ways. The principal improvement has been in the 
ability to use more grid points the original 64 x 64 x 64 grid and, in some 
recent calculations 128 grid points have been used in one or two of the direc- 
tions. They have also made improveifients in the subgrid scale model and in the 
numerical method. 

Kira and Moin (1979) reported the results of 64 x 64 x 64 simulations 
with a model which damped the subgrid scale viscosity near the wall more 
strongly than the previous model. We shall look at some of their results. 
The mean velocity profile they obtained is compared with several experiments 
in Fig. 7.3. 'Fhe existence of a logarithmic region in the flow with the 
correct slope is one of the major achievements of the whole of higher-level 
simulations. The profile near the wall lies below the expected profile (indi- 
cated by u = y in the figure), and this is probably due to a weakness of 
the model in the region near the wall. The components of the turbulence are 
shown in Fig. 7.4, although the experimental data are not shown, the agree- 
ment is quite good — the experimental data show quite a bit of scatter. The 
pressure-strain terras are shown in Fig. 7.5. In the center o£ the channel. 
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these terms drain energy from the fluctuations of the streaiuwise velocity and 
transfer it to the other components; this accords with expectation. However, 
near the wall there is a large transfer from the normal fluctuations to the 
spanwise ones. This had also been found by Moin et al. (I97U) and was unex- 
pected. It is apparently due to fluid moving toward the wall being stopped by 
the wall. The vertical motions are converted into horizontal motions, and the 
result of this "splat" effect and the normal energy transfers is shown in the 
figure. Again, the quantitative results may be incorrect, but it is unlikely 
that the qualitative result is incorrect. More recent (and more accurate) 
results by Moin and Kim (19B1) show a smaller, but still significant, "splat" 
effect. 

Contours of the fluctuating velocity on a plane parallel to and close to 
the wall are shown in Fig. 7.6, the presence of long streaks is obvious. A 
similar plot for a plane close to the center of the channel is shown in Fig. 
7.7; there is no evidonce of streaky behavior at this plane. A number of 
other plots of this kind are given in their paper. 

In a more recent paper, Kim and Moin (1981) have done calculations with 
still greater resolution and further improvements in both models and numerical 
methods. The results are qualitatively similar to those presented above but 
differ quantitatively. Hiey have also produced a simulated flow-visualization 
motion picture that duplicates most of the phenomena observed in laboratory 
motion pictures. Tliis application of the results should play a very important 
role in the future. 

The splat effect is also seen in the shear-free wall layer. This is 
simply a turbulence near a wall whicVj is moving at the same mean velocity as 
the wall. The precise nature of this flow depends on the Reynolds number. A 
simulation by Biringen and Reynolds (1981) captured most of the effects ob- 
served in the experiments. However, we shall not review these results in 
detail here. 

3. Transition 

As stated in the introduction, transition in boundary layers is a subject 
of great technological importance. However, transition is very sensitive to a 
number of factors, including the precise velocity profile, the level of the 
disturbance, wall roughness, etc. As a result, transition experiments are 


83 


very difficult to pt!>rform reproduceably , and there is considerable scatter in 
the data. Naturally, simulations of these flows will be very sensitive to 
similar factors. Hius , a great deal of care will be necessary to simulate 
these flows. 

Ldnear stability theory predicts that the laminar boundary i ayer profile 
is unstable with respect to disturbances that result in Tollmieii-Schlichting 
waves, Tliis instability is much 'ss explosive than tliat of the free shear 
layer. It is generally believed that the Tollmien-Schlichting waves grow 
until nonlinear effects take over and complex interactions lead to the fully 
turbulent boundary layer. However, the late stages of transition are poorly 
understood. 

The first direct simulation of transition in wall-bounded flows was made 
by Kells and Orszag (1979) and Orszag and Patera (1980,1981). They chose to 
study channel flow at Reynolds numbers for which the flow is linearly stable. 
However, transition does take place at the Reynolds numbers studied. In their 
simulation, Orszag and Patera took a fully developed laminar channel profile 
(Poiseuille profile) and added finite-amplitude two-dimensional Tollmien- 
Schlichting waves to it, these viaves are different in the channel than in the 
boundary layer. Tliey found that the waves decayed slowly and that the rate of 
decay decreases as the Reynolds number increases, this is expected. However, 
they found that, when a small three-dimensional disturbance is introduced into 
the flow, it grows very rapidly. TSie growth of the three-dimensional wave is 
rapid enough to enter the nonlinear regime before the two-dimensional wave has 
decayed. At this point the simulation develops considerable energy at high 
wavenumbers and has to be stopped; as there is no model in the simulation, 
there is no way to continue. However, this simulation has provided an expla- 
nation of the instability of this flow; it is apparently due to the three- 
dimensional instability of stable two-dimensional waves. Orszag and Patera 
(1981) have done similar simulations for Couette and cylindrical tube flows. 

Some of their results are shown in Fig. 7.8. Hie decay of the two- 
dimensional wave and the growth of the three-dimensional wave are quite 
apparent. 

A simulation of the instability of the boundary layer has been made by 
Wray (unpublished). In order to avoid the difficulty that arises from the 
spatially developing boundary layer, he chose a time-developing boundary 
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layer; physically, this corresponds to the boundary layer that develops on a 
suddenly started plate. Although the velocity profile of the time-developing 
boundary layer is different from that of the spatial layer, the calculation 
was started with the Blaslus profile appropriate to the spatial layer. To 
this profile, a weak Tollmien-Schlichting wave and a weak three-dimensional 
random disturbance was added. 

The disturbance grows very slowly at first (as expected) until it builds 
up to a level at which nonlinear effect!^ become Important. At this point, the 
rate of change of the layer becomes spectacular. Tlie contours of the various 
velocity components and the vorticity develop more and more structure. Compa- 
risons with experimental results for the parameters reveal a considerable sim- 
ilarity; the comparison is necessarily qualitative, but it is remarkably good. 

Eventually, Ichis simulation developed a considerable amount of energy at 
high wavenumbers, and it had to be stopped. There is no way to continue this 
simulation beyond this point without more resolution. Unfortunately, it may 
be that the small scales play an Important role in the development of this 
flow, and it is not known whether the addition of a model will cure the prob- 
lem. 
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Chapter VIII 
OTHER APPLICATIONS 


As we have stated earlier, higher-level simulation began in meteorology 
and oceanography. These fields have maintained an active Interest in the 
simulation of the global circulation of the Earth's atmosphere and oceans. 
TIte methods used are similar to the ones described in this report, but there 
are additional difficulties. The principal of these is that thermal energy 
and the transport of water vapor (in the atmosphere) and salt (in the ocean) 
are very Important in these flows, and one must deal with the effects of 
stratification, evaporation, and condensation. When coupled with the limita- 
tion to very coarse grids '200 km Is typical in these Simulations today) , we 
see that the problems are considerably more difficult than the ones dealt with 
in this report. They are, however, of great importance, and considerable 
effort is being devoted to them. Die author has only a passing knowledge of 
the Work in these areas, and this is the reason why the subject is not covered 
in this report. 

Methods similar to the ones given in this report have also been applied 
to smaller-scale environmental problems. For example, simulations have been 
made of local parts of the atmosphere by these methods; these are called 
mesoscale simulations. The author is familiar only with a few papers by 
Iteardorff in this area, in these papers, he used a complete Reynolds stress 
model for the subgrid scale Reynolds stresses. Others have applied these 
methods to the prediction of the flow in lakes, harbors, and other small 
bodies of water. Of the work in this field, the author is familiar only with 
some of what has been done at his institution. Findikakis (1960) has recently 
developed a finite-element method for computing such flows. 

On a still smaller scale, there, have been a number cf extensions of the 
work covered in the earlier sections of this report. Schumann and his cowork- 
ers have used the method that was described for channel flow for the simula- 
tion of flows in annuli and made other extensions. In particular, they have 
computed the channel and the annulus with heat transfer, in the simulations, 
the temperature is treated as a passive scalar. We have not dealt with this 
Work at length in this paper for several reasons. It is covered in detail in 
the report cf Schumann et al. (1980). Also, since the results produced by 
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Schumann’s method differ consldecabiy from those of Kim and Motn (1979) and 
Moin and Kim (1981) for the channel flow, we are unsure about the accuracy of 
the method when applied to heat transfer. For similar reasons we have not 
covered their work on the effect of roughness. 

Schumann, Crotzbach, and Klelser have applied their method to natural 
convection flow between parallel horizontal plates. They covered a very large 
range of Rayleigh number and were able to predict the observed transitions 
from one flow regime to another. lliis is an excellent piece of work and was 
not covered because it did not fit any of the subject tieadings used in tills 
paper. Crotzbach (1979) has also Investigated simulated flows in vertical 
channels with the influence of buoyancy. 

Finally, we shall mention a method that competes with the, grid-based 
methods that are the primary subject of tills report. 'Hiese are methods in 
which the motions of vortices are followed (vortex-tracking methods). A number 
of interesting features of transitional and turbulent flows have been computed 
by this method, including flows with separation. The full capabilities of 
this approach and a comparison of it with the methods discussed in this report 
are given in a review paper by Leonard (1981). Hybrid methods which use some 
ideas from vortex-tracking and some from grid-based methods are also being 
Investigated at the present time, the interested reader is referred to the 
paper by Couet, Leonard, and Buneman (1980). 

Tliere ire undoubtedly areas that have been overlooked in this report. 
The author has tried his best to be complete, but in any subject area that has 
become this large something is likely to be missed. Tliere is no intent to 
minimize any contributions that have been missed. 


87 


Chapter IX 

CONCLUSIONS AND FUTURE I'iRECTlONS 


1 . Where Are We Now? 

We hope that this report has shown that higher-level simulations of tur- 
bulent flow have reached a point in their development which allows them to 
play an important role in turbulent fluid mechanics. Let us now sum up where 
the field stands today. We start with the positive points. 

a) The basic ideas of large eddy simulation seem sound. Specifically, 
they seem to be able to handle homogeneous turbulent flows and free shear 
flows quite well. For wall-bounded flows, the importance of small structures 
near the wall is a problem, and these flows are difficult to deal with, but 
gcod progress has been made. 

b) Direct simulation of many interesting flows are now feasible. We are 
limited to low Reynolds numbers, but this restriction may not be important in 
some flows, as the large scales may be nearly Reynolds number independent. 
Alternatively, one can regard the viscosity as a simple subgrid scale model 
and pretend that a higher Reynolds number flow is being simulated. Both of 
these approaches have been taken. Orszag has used the concept of "Reynolds 
number similarity" with considerable success. Rubesin (1979) regarded direct 
simulations as large eddy simulations, also with ccasiderabie success. 

c) Higher-level simulations have come to the point at which they are 
able to provide information on quantities that are difficult to measure in the 
laboratory. In this role, they are able to evaluate turbulence models in a 
way that is very difficult to do by any other method. 

d) Higher-level simulations are able, in some cases, to simulate flows 
that are very difficult or impossible to perfom in the laboratory. Some 
examples are flows with rotation and/or compressibility. 

e) It is now possible to do flow visualizations using full and/or large 
eddy simulations. Hiese visualizations reproduce much of what is seen in the 
laboratory. Tliey also offer flexibility that is difficult to match in the 
laboratory- They can be used to look in detail at specific regions, and can 
even be used backwards in time. 
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Now let 46 consider some of the difficulties. 

a) The most obvious problem is that these methods require large amounts 
of computer time. Although some of the simpler flows can be done in a few 
minutes on large machines, running times of the order of hours arc not unusual 
for the more difficult flows. Use of these methods must be restricted to 
individuals with access to the machines that can do these simulations. Some 
means of assuring that the problems of greatest interest are done is 
necessary. 

b) Although some flows are amenable to full simulation, Reynolds number 
similarity does not hold for all flows, so it is not possible to treat low 
Reynolds number flows as models of a high Reynolds number flow in all cases. 
Better subgrid scale models will be necessary if high Reynolds number flows 
are to be simulated, but it may be very difficult to find models with wide 
applicability . 

On balance, the contribution of higher level simulations seems to be more 
than worth the cost, and the approach is just beginning to reach its poten- 
tial. With new generations of computers, it should be possible to do much 
more with these methods. 

^ • Where Are We Going? 

It is clear that a great deal remains to be done in turbulence simula- 
tion. 'fliere are many directions which can be taken in the future, and, with 
more group, s beginning to do these types of simulations , we expect the area to 
broaden rapidly. Of course, it is difficult to predict the future with any 
precision, but it is always interesting to try. Let us look at what can be 
expected in the next few years. 

a) One obvious direction in which highei-level simulations will be 
extended is toward the simulation of a larger number and greater variety of 
flows. There are many possibilities, so the following list cannot be all- 
inclusive. 

i) The flows which have already been simulated can be done with 
additional effects. Hius , to any of the flows treated in this report we can 
add rotation, curvature, heat transfer, passive scalars, and/or pressure 
gradients alone or in combination. In the wall-bounded flows we can also add 
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wall roughnasB and blowing oc succlon* Many of tlteae cffectB are quite 
Impoctv'tnt in engineering flows and should be considered at an early date. 

11) To date, no method has been found for dealing with Inflow or 
outflow boundaries. Tlie outflow boundary can probably be handled by the usual 
method of requiring the streamwlse derivatives to be zero at the outlet. TYie 
inflow condition is much more difficult, because it is necessary to prescribe 
a realistic representation of the turbulence in order not to require too much 
of the computation; to do this would waste a very large pact of the computa- 
tional resource. Being able to handle inflow and outflow boundaries is cen- 
tral to the computation of many flows of interest. 

ill) There are some fairly simple flows which have not been done. 
Included among these are the jet and the wall jet. 

b) Simulation of wall-bounded flows is much simpler if the kinds of 
boundary conditionc used by Deardorff and Schumann can be applied. Chapman 
(1980) estimated that the savings to be realised In this way could make the 
difference between practical use of the higher-level simulations and their 
continuing to be confined to research. Accurate boundary conditions of that 
type need to be searched for. 

c) Use of higher-level simulations in conjunction with flow visualiza- 
tion and statistical methods should become a very powerful tool for investi- 
gating the structure of turbulent flows. It is possible that such an approach 
may be able to tie the structure of turbulent flows to the modeling. This is 
highly speculative, but, if it can be done, it could be an important step for- 
ward. We may become "computational experimentalists." 

d) The interaction of higher-level simulations and conventional model.- 
should become stronger. We can envision a time when people developing new 
models will routinely validate them using higher-level simulations. Cer- 
tainly, we can expect higher-level simulations to be helpful in determining 
the constants in the models. It is worthwhile to set up a facility which is 
available for this purpose. 

e) We expect that there will be considerable work on the improvement of 
subgrid scale models, but the direction this work will take is not obvious. 

f) Higher-level simulations will be extended to include a number of phe- 
nomena chat are not currently treated. Sound generation should be relatively 
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eafiy , as It seems to depend mainly on the large 8csles< Combustion should be 
very challenging, because the chemical reaction depends on intimate mixing at 
the small scales • 

g) Something has to be left to the reader's imagination. 
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Poutler Space 



a) Box Filter 



b) Fourier Space Sharp- cutoff 




c) Gaussian 


Figure 2.1 Some Filters Commonly Used in Large Eddy Simulation 
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Control volume £or x-momcvitum equation 





h — Control volume for y- 
I momentum equation 
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Figure 2.2 The Staggered Crld and the Control Volumes. 
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(a) 


(b) 


Figure 2.3 Definition of the Ifcge scale field by 

(a) the Deardorff averaging method 

(b) filtering 
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EXACT VALUE 
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sr.'n 


“^00 -250 ^200 -IS)U -iw - 

MODEL VALUE 


Figure 3.1 Scatter plot of th7s!Irotrnir">^dIl'°in 

Is 0.4 Cor this case. From mcnxj.Aa 




Figure 3.2 


The constant in the Smagorinsky nodel as a function of the sungrid 
scale Reynolds number R^^g * IsfA^/v where jsj is the r.m.s. 
large scale strain. From McMillan and Ferziger (1978) . 



k 


Figure 3.3 A typical turbulence spectrum and the effect of filtering. 

Region 1 represents the smallest scales of the filtered 
or resolved component of the field. Region 2 represents 
the largest scales of the subgrid scale or unresolved com- 
ponent of the field. The scale similarity model relates 
these. 
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Figure 3.^4 Scatter plot of u^3Tj^j/3Xj (the dlsHlpatlon due to 

the aubgrid scale model) for the scale’-slmllarity model 
in strained turbulence. The correlation is Improved 
relative, to the model shown In Figure 3.1 but the. model 
has no average dissipation. From McMillan et al, (I960). 
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Fluure A.Ih The effective wave membei; for a number of approximations. 
The methods used are; kA , pseudospectral ; k"A , 
order central difference; k’A , fourth order central dif 
ference; kA , compact method, 
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Figure 4.1b The square of the effective wave number in various 
second derivative approximations. The methods ate 
given in Figure 3.1. From Mansour et al. (1978). 
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Figure 5.1 The turbulent kinetic energy and its 
components in a direct simulation of 
a turbulent flow at r 40. 


Figure 5.2 The skewness history in the flow of 
Figure 3.1. 


From Feiereisen et al. (1981) 
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Figure 5.4 The decay Index for isotropic turbulence as a function 

of the Reynolds number. Crosses are computational results; 
other symbols are experiments. From Shirani et al. (1981). 
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Figure 5.5 The skewness of the velocity derivative as a function of Reynolds nuaber for Isotropic turbulence 
Symbols are: x direct simulations; + large eddy simulations; all other symbols are experi- 
mental data. From Shiranl et al. (1981). 
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Figure 5.6 



The filtered energy spectrum from a typical large eddy 
simulation. Points are computed results; the .,ine 
represents the filtered experimental data. 
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Figure 



5.7 A number of quantities on a direct simulation of anisotropic 


turbulence. 

The quantities are: 


Eij : 

Reynolds stress tensor 


eij : 

Dissipation tensor 


4>ij : 

Pressure strain tensor 


^ij = 

Integral length scales 



The fifth is the r.m.s. pressure gradient 
sixth, the r.m.s. pressure fluctuations. 

Schumann and Herting (1976) . 

and 

From 
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Figure 5.8 Rotta's constants C and C' (see eqs. 5-6, 5-7) 
versus time for several axlsynunetric anisotropic 
cases. From Schumann and Herring (1976). 






PECAV Of TURBULENCE 



Figure 5 
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.9 The decay of rotating turbulence. Lines are 
large eddy simulation; points are data. 
Rotation was present from the 
calculation, from Bardina et al, (19B1;. 
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FiKure 5.10 Direct simulation of rotating turbulence. The 
flow was all.<wed to relax before rotation was 
"turned on." From Bardina et al. (1981). 
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/3t/^ 


Figure 5.11 The behavior of turbulence undergoing 
strain. From McMillan et al, ( 1980 ) . 
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St * 0.5 


St " “0.5 


Figure 5,12 The mesh used in sheared turbulence. Flow is 
started with coordinate system shown on top. 
When it has sheared to the position of the 
middle figure, the flow is Interpolated on to 
the grid shown at the bottom. From Feiereisen 
et al. (1981) . 
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Figure 5.14 The energy spectrum and its components at two 
different times in sheared turbulence. From 
Feier eisen et al . (1981) . 
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Figure 5.16 The spectra of the various terms that contribute 
to the change of the 3-D energy spectrum. From 
Shirani et al. (1981). 
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Figure 5.17 The dissipation anisotropy vs the Reynolds stress 
anisotropy in sheared turbulence. From Peiereisen 
et al, ( 1981 ). 
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Figure 5.18 


The Rotta "constant" C in sheared 
turbulence. The abscissa is the 
model constant obtained from the simu- 
lations. The ordinate is the ratio of 
the abscissa to the value obtained from 
a least squares fit. A perfect fit 
would give a value of 1.0 for the ordi- 
nate; the values for the various compo- 
nents have been shifted for clarity. 
From Feiereisen et al. (1981) . 
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Figure 6.1 A schematic of the mixing layer. 
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Figure 6.2 Momentum thickness of a developing mixing layer 
vs time; low initial intensity cases. From 
Cain et al. (1981). 
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Figure 6.3 Momentum thickness of a mixing layer vs time; medium | 

initial intensity cases. From Cain et al. (1981). 
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Figure 6.4 Momentum thickness of a mixing layer vs time; high initial 
Intensity cases. From Cain et al. (1981). 
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Figure 6.5 Turbulence intensity at center of mixing layer vs time; 
low Intensity cases. From Cain et al. (1981). 
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Figure 6.7 


Turbulence intensity at center 
high initial Intensity cases. 


of mixing layer vs time; 
From Cain et al. (1981). 
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Flaure 6.10 "Dye lines" late in mixing layer development; 
* low Initial intensity case. From Cain et al. 

/ 1 Qfti \ . 
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T=171 


Figure 6.11 


"Dye lines" late in mixing layer development; medium 
initial Intensity case (same as Figure 6.10 except 
for intensity). From Cain et al. (1981). 
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Figure 6.12 Decay of Maximum Mean Velocity (Uj,) and Maximum Axial 
Turbulent Intensity (u„) in From Riley and 

Metcalfe (1978). 
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Figure 6.13 


Growth of Mean 
Radius (rx). 


Wake Radius (r,„) and Turbulent Wake 
From Riley and Metcalfe (1978). 
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Figure 7.2 The 'constant' in the model of the pressure-strain term 
as a function of the normal coordinate. From Schumann 
et al. (1981). 
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Figure 7 . 3 The 
of I 


Figure 



.4 Turbulence intensities In channel flow (resolved component 
only) . From Kim and Koln (1979) . 
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7.6 Contours of fluctuations of the streawise velocity 
near the wall. Froa Kia and Mo in (1979). 
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Figure 7.7 Contours of the fluctuations of the streaawlse velocity 

near the center of the channel . From Kim and Moln (1979) . 
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Ener«y of t«- .od three-dli.enolon.1 «aveo 

UmlMt channel flow. The 2-D >«»« 

From Orszag and Patera (1981) . 
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